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ABSTRACT 

Although it is known that the stalled accretion shock in models of core-collapse supernovae turns 
into explosion when the neutrino luminosity from the proto-neutron star (PNS) exceeds a critical value 
core) (t ne "neutrino mechanism"), the physics of L^ rl * ore has never been systematically explored. 
We solve the accretion problem between the proto-neutron star (PNS) surface and the accretion shock. 
We quantify the deep connection between the general problem of accretion flows with bounding shocks 
and the neutrino mechanism. In particular, we show that there is a maximum, critical sound speed 
above which the shock jump conditions cannot be satisfied and steady-state accretion is impossible. 
This physics is general and does not depend on a specific heating mechanism. For the simple model of 
pressure-less free-fall onto a shock bounding an isothermal accretion flow, we show that shock solutions 
are possible only for sound speed ct < c" 1 * and that c^/v^ sc — 3/16 = 0.1875 at c^ 1 *. We generalize 
this result to the supernova problem, showing that the same physics determines L£, rl * oro . The critical 
condition for explosion can be written as c|/Wg SC ~ 0.19, where eg is the adiabatic sound speed. This 
"antesonic" condition describes L"^ OIC over a broad range of parameters, and other criteria proposed 
in the literature fail to capture this physics. We show that the accretion luminosity reduces L"^ m . c 
non-trivially. A larger PNS radius decreases £™* ore , implying that a stiff high-density equation of 
state may be preferred. Finally, using an analytic model, we provide evidence that the reduction 
of L"c 0re seen in recent multi-dimensional simulations results from reduced cooling efficiency, rather 
than an increased heating rate. 

Subject headings: Accretion, accretion disks 
supernovae: general 
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1. INTRODUCTION 

When the core of a massive star collapses at the end of 
its life as a result of the Chandrasekhar instability, the 
central matter reaches nuclear densities, stiffens dramat- 
ically as a result of the hard core repulsion of the strong 
force, and drives a Shockwave into the super-sonically col- 
lapsing outer mantle. The outward progress of the shock- 
wave is halted by a combination of energy losses from the 
dissociation of nuclei across the shock, neutrino emission 
as the shockwave moves from the optically-thick inner 
core to the optically-thin surrounding envelope, and the 
ram pressure of the overlying inf ailing outer iron core. 
A quasi-static accretion phase then ensues, which lingers 
for many dynamical times before eventual explosion. It 
consists of a hot, optically-thick, accreting proto-neutron 
star (PNS) with neutrinosphere radius ~ 30 to 60 km, 
radiating its gravitational binding energy in neutrinos of 
all flavors, surrounded by a standoff accretion shock with 
radius of ~ 200 km. 

Despite decades of modelling effort, the mechanism 
responsible for reviving the shockwave to positive 
velocities has not been conclusively identified. The 
"delayed neutrino mechanism" relies on neutrino 
heating in the semi-transparent subsonic accretion 
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flow to drive explosion. A fraction of the neutri- 
nos diffusing out of the collapsed PNS deposit their 
energy in the accretion flow, and for a sufficiently 
large heating rate explosion results . This mecha- 
nism was originally discusse d bv iColeate fe White! 
(119661 ). and then developed bv lBethe fe Wilson! (|1985T T 
However, the most detailed one-dimensional mod- 
els generically fai l to explode, except f or the lowest 
mass p rogenitors (iRampp fe Janka! 120001: iBruenn et al. | 

IMezzacappa et all 
Kitaura et al.l 120061: 



2001 



2001 



Janka et al.l 



ILiebendorfer et al l 1 2001 
iThompson et al.l I200l ~ 



120081 ). These models clearly 
multi-dimensional effects such as convection 



neglect 
both in 

the PNS interior (Keil et al.l 119961 : IMezzacappa et ah! 
[19981 iDessart et al.l I2006D . and in the h eating region 



betw ee n the ne utrinosphere and t he sho ck ( Hera nt et al 
1991 Uanka fe Muller! 119961: IBurrows et al 



Frver fe Hegerl l2000t iFrver fe Warred I2002L 



shock 
instability 



frver (z ricge n izuui 

iBuras et al.l l2006aft . the standing accretion 
insta bility (S ASI) and/or vortical- acoustic 



( Fogl izzd 



Blondin fe Mezzaca ppa 



Blondin et al.| 


2003; 


Ohnishi ct al. 


2006; 



Iwak ami et al.1 
2010f ). and the potential for other multi-dimensional 



n on-linear phenomena such as the "acoustic" mechanism 
of IBurrows et al.| ( [2006. 200 71), which was critically eval- 
uated bv I Weinberg fe Quataertl (|2008l ). Indeed, many of 
these works have claimed that multi-dimensional effects 
are crucial for explosion because convection between 
the PNS and the shock allows the matter to stay in 
the region of net energy deposition longer, and an 
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increase in the average shock radius as a result of the 
SASI and convection makes the heating region larger, 
while putting more matter higher in the gravitational 
potential well. 

An approach complimentary to these multi- 
dimensional modelling efforts was followed by 
iBurrows fe Goshvl (|1993l ) who calculated the steady- 
state structure of the region between the PNS and the 
accretion shock, during the quasi-time-steady epoch 
following collapse, bounce, and shock formation. For a 
fixed mass accretion rate, they increased the PNS core 
neutrino luminosity by hand until no steady-state solu- 
tion was possible. The critical core neutrino luminosity 
(i"core) that bounds steady-state accretion solutions 
from above has been identified with the initiation of the 
supernova explosion. With this knowledge in hand, the 
task of getting a supernova simulation to explode via 
the neutrino mechanism can be understood as either 
an effort to increase the neutrino flux (or effective 
heating rate) up to the critical value, or to decrease 
^"corc so tnat can b e reached with the available 
lu minosities. Examples o f the former include the work 
bv lWilson fc Mavlel (|1988| ). which appealed to salt-finger 
(doubly-diffusive) convection within the PNS to increase 
the core lumino s ity a t early times after collapse, or 
I Thompson et al.l ((20 05) , who found that viscous heating 
in models of rotating collapse co uld enable explosion. 
Exam ples of the latter include Ya masaki fe Yamadal 
(2006), who showed that convection in the heating re- 
gion generically decreases £ "'* ore , or lMurphv fc Burrows! 
(|2008fl and iNordhaus et al.l fcoiOt ) , who found the sa me 
in 2D and 3D simulations (see also lHanke et al.ll2011l ). 

However, it is to a great extent unknown why there is 
a critical neutrino luminosity at all, what determines its 
existence, and how it scales with the parameters of the 
problem. Indeed, the question of the physics of £™c Or0 
is particularly paradoxical because during the stalled ac- 
cretion shock phase, the size of the PNS, its neutrino lu- 
minosity, and the mass accretion rate through the shock 
change significantly, yet the simulations show that the 
position of the shock moves only very slowly. This im- 
plies there are many possible near-equilibrium configu- 
rations that yield similar shock radii for very different 
global parameters. For example, in low mass progeni- 
tors, the incoming mass accretion rate can decrease by a 
factor of 10 in the first second after bounce, but the neu- 
trino luminosities vary by less than a factor of 2, and yet 
the shock radius is essentially constant. Given this early- 
time behavior, if the neutrino luminosity were to increase 
very slowly at a later time in the evolution, why would 
one expect a catastrophic change from steady accretion 
to dynamical outward expansion? One might instead ex- 
pect that an increase in the net energy deposition would 
simply result in the shock just getting pushed to some 
different equilibrium position with larger shock radius. 
Contrary to this expectation, the existence of L"^ OIC 
implies instead that if the energy deposition below the 
shock reaches a critical value, the shock can no longer ex- 
ist together with steady-state accretion, the Shockwave 
expands dy namically, and the star e xplodes as a super- 
nova (e.g., IBurrows fe Goshvl 119931: iHerant et "all 119941: 
IBurrows et al.lll995l : Uanka fe Mullen 1 1 1 9961 ) . 

Within the framework of time-steady spherically sym- 



metric models of supernovae, this paper provides an 
answer to two questions connected to these problems, 
and a possible explanation for a third. First, what 
is the physics of L™* orc ? In Section [2j we investi- 
gate isothermal and polytropic accretion flows with a 
bounding Shockwave and we show that even in such 
a simplistic setting there exists a critical value of the 
controlling parameter — the isothermal sound speed 
- that separates steady-state accretion solutions from 
wind solutions , which are identified with explosions 
(IBurrows! 119871: IB urrows fe Goshvl [1991 IBurrows et all 
119951 : lYamasaki fc Yamadal I2005T L This critical value 
is determined by conservation of momentum and en- 
er gy across the standoff accretion s hock as was found 
by lYam asaki & Yamada] (|2005l 12006!): above the critical 
value it is not possible to satisfy these conditions with 
a steady-state accretion solution, and instead the flow 
must rearrange itself into an outgoing wind. In Section[3l 
we calculate the steady-state structure of the accretion 
flow between the neutrinosphere and the shock in the 
supernova problem. We explicitly s how that the criti- 
cal ne utrino luminosity discovered bv IBurrows fc Goshvl 
(1993) is equivalent to the critical condition from the 
isothermal model. That is, the neutrino mechanism of 
supernovae as formulated bu \Burrows & GosM IBM) 
identical to the statement that the shock jump conditions 
cannot be satisfied together with the Euler equations for 
steady-state accretion when the sound speed of the flow 
exceeds the critical value. Furthermore, we discuss the 
effects of radiation transport and the accretion luminos- 



ity on the absolute value of L 



crit 



and we show how the 



latter depends on the mass and radius of the PNS and 
the mass accretion rate. 

The second question this paper answers is "What cri- 
terion does £™' ore correspond to in terms of variables 
within the accretion flow itself?" For example, it has 
been claimed that when the advection time in the heating 
region becomes l onger than the heating timescale, explo- 
sion must result ([Ja nka 2001 : [Thompson fc Murravll2001l: 
Thompson et al.l 120051: iBuras et al.l l2006bt iScheck et al.l 
20081: iMurphv fc Burrows! 120081 ). Is this heuristic crite- 
rion identical to the critical neutrino luminosity? The 
answer is no. In fact, we show in Section |4] that L"' l ^ OIC 
corresponds to a virtually constant ratio of the sound 
speed to the escape velocity in the accretion flow — the 
"antesoni<0 condition:" c| ~ 0.l9v^ sc , where v csc is the 
escape speed — and that this condition is a consequence 
of the physics of the critical luminosity itself. 

Third, an intrig ui ng observation was mad e by 
l Ohnishi et all (120061) IMurphv fc Burrows! (|200l and 
INordhaus et al.l (|2010l ) who found that the mere fact of 
increasing the dimension of the simulation lowers L£, rl * ore . 
While going from ID to 2D allows for the effects of con- 
vection and the SASI, it is not clear what is gained by 
going from 2D to 3D, or why L^ OTe should decrease 
monotonically with dimension. In Section [5l we provide 



evidence that the reduction of if 



observed in 2D 



4 We call this condition "antesonic" because the required sound 
speed is significantly below the local escape velocity, and the critical 
radius at which this point occurs is smaller than the radius of sonic 
point in Bondi accretion. This condition is reached in supernovae 
at the time of explosion, before the flow arranges itself into a super- 
sonic neutrino-driven wind (see Section [2] for details). 
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and 3D simulations occurs because the flow can orga- 
nize itself to cool less efficiently and yet still maintain 
hydrostatic equilibrium, giving an overall lower cooling 
efficiency, lower critical luminosities, larger entropy, and 
larger shock radii. Thus, the reduction of L c J l l me m 2D 
and 3D may be due to less efficient cooling, rather than 
more efficient heating. 

2. ISOTHERMAL ACCRETION BOUNDED BY A SHOCK 

Here we review the basic physics of spherically- 
symmetric isothermal accretion with a shock. We find 
that this idealized model problem is a key to under- 
standing the mechanism of supernovae. We also discuss 
here the analogous problem of polytropic accretion flows 
with a constant adiabatic index T before solving the more 
complete supernova problem in Section 3. 

2.1. Topology of solutions 

The velocity structure of spherical steady-state isother- 
mal flows is described by the equation 

V Jf ) dx x 2x 2 ' w 

where ^# = v/ct is the Mach number, v is the fluid ve- 
locity, ct is the isothermal sound speed, x = rc T /(2GM) 
is the rescaled radial coordinate r, and M is the mass of 
the central object. It is possible for a standing shock 
wave in the flow to exist at a point that satisfies the two 
Rankinc-Hugoniot shock jump conditions 

p-j#-=p + Jt+, (2) 
p-{JK-) 2 +p-=p + (^ + ) 2 +p+, (3) 

which express conservation of mass and momentum, re- 
spectively. Here, p is the mass density and the + and 
— superscripts correspond to the quantities evaluated 
just upstream and downstream of the shock, respectively. 
The physically relevant solution to equations ([2]-[3]) is 
.//'.// = 1. 

In Figure [TJ we show with black solid lines Mach- 
number profiles of flows with a single constant mass ac- 
cretion rate M = — 1 M Q s _1 obtained by solving equa- 
tion ([T]). The integration starts at a fixed radius r v cor- 
responding to the surface of the star, and with a fixed 
velocity v v — M /(Aicr^pv), where p v is the prescribed 
mass density at r„. Different lines correspond to differ- 
ent values of c T . The dashed line shows the decelerating 
transonic solution, which goes through the sonic poinlQ 
located at atonic = 0.25 and ^# son ic = — 1< It has high 
velocity at large radii and j% — > as x — > 0, unlike Bondi 
accretion flow, which starts with zero velocity at infin- 
ity, goes through the sonic point and has ^# — > — oo as 
x —> (shown with dash-dotted line). The decelerating 
transonic solution separates accretion "breezes" , which 
are always subsonic (— 1 < jM < 0) and denoted as B, 
from solutions going through j$ = — 1 and denoted as 

5 In this paper we assign the name "sonic point" to the position 
of the flow where the numerator and denominator of the fluid mo- 
mentum equation simultaneously vanish. Although this point is 
generally called the "critical point" , we avoid this name to prevent 
confusion with the critical neutrino luminosity discussed later in 
the paper. 



A. Although the A-type solutions lying below the decel- 
erating transonic solution are considered unphysical in 
standard Bondi accretion theory, the parts of these so- 
lutions that span from r„ to the shock are viable in this 
setting (|McCrealll956h . 

The dotted lines in Figure Q] show the velocities just 
downstream of the shock, corresponding to three 

different assumed profiles of the upstream velocities ^ + : 
(1) Bondi flow (blue), (2) free fall with Jl+ = -x" 1 / 2 
(green), and (3) pressure- less free fall (red), that is a good 
approximation in more realistic supernova calculations 
([Colgate &: White! I1966H and that requires modification 
of equation ([3]) by assuming that = — a; -1 / 2 <C — 1, 
and 

p-(.^-) 2 + P - = p+ x - 1 . 

The physically relevant solution to equations © and (|3"aT) 
is J(~ = {\Jx~ x — 4 — x~ 1 / 2 )/2. For the sake of clarity, 
we do not show the upstream profiles in Figure [1] except 
for Bondi flow (dash-dotted line). 

The velocity profile of the shocked accretion flow for a 
given ct is constructed by following one of the black lines 
from r v until it crosses a dotted line that corresponds to 
the appropriate shock jump conditions. At this intersec- 
tion, the solution jumps to whatever velocity profile 
was assumed to be present upstream. 

The essential point of FigureQ]is that for a fixed M the 
presence of a shock is not guaranteed for all values of ct- 
For example, consider the blue dotted line, which corre- 
sponds to Bondi accretion flow upstream of the shock. 
For the smallest values of Ct, the required shock radius 
— the intersection of the black solid line with blue dotted 
line — would be below the radius of the star, an unphysi- 
cal situation. As we step to higher Cy, the shock appears 
at r v and moves outward. At a critical value c™ cor- 
responding to the decelerating transonic accretion flow 
(dashed line), the shock can only coincide with the sonic 
point at x mt = atonic = 0.25. Here there is no jump 
in the velocity, and the shoc k degenerates into a jump in 
the derivative of the velocity (Velli 1994; Del Zanna et al. 
1998). A shock is not possible for ct > c^ rlt , because the 
blue dotted line lies below the decelerating transonic so- 
lution (dashed line) in the range of radii of interest and 
hence B-type solutions are not viable solutions inside the 
shock - that is, solutions of type B do not intersect the 
blue dotted line. 

The situation is somewhat different when the flow up- 
stream of the shock is in free fall, as shown with the red 
and green dotted lines in Figure [TJ In these cases, the 
shock jump conditions allow part of the B-type solutions 
as viable solutions downstream of the shock. Further- 
more, for values of ct higher than that corresponding 
to the decelerating transonic solution (dashed line) two 
shock radii are possible for a single value of ct, because 
the curvature of B-type solutions and shock jump condi- 
tions yields two intersections. When stepping to higher 
values of ct, the two shock radii come closer together and 
finally coincide at a critical value c™'. Contrary to the 
case of Bondi flow above the shock, there is no relation 
between the critical sound speed c™' and the sonic point. 
The shock radius at c" 1 * also does not coincide with the 
velocity minimum of any of the B-type solutions. Above 
the critical sound speed c" 1 *, a shock in the flow is not 
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Figure 1. Isothermal accretion plotted in the space of Mach number and rescaled radial coordinated x = rc^,/(2GM). Solid black 
lines show solutions to eq. JT} with M = —IMq s~ 1 and M = 1.4 Mq starting from rv = 30km (grey dashed line), and with fixed 
velocity v„ = M /(Airr^p^), where p v = 3 X 10 10 g cm 3 . The value of eL increases from 4 X 10 18 (black solid line starting at lowest x) to 

D 17 cm 2 s — 2 . The dashed line is the decelerating transonic solution going through 
0.25. The dash-dot line is Bondi accretion flow. Dotted lines are velocities just downstream of a 



1.' 



x 10 1 



" ]_ clt Xa, 



the sonic point ~4 

shock positioned at any x, assuming that the upstream flow is either Bondi accretion flow (blue bottom line), in free fall (green top line), 
or in pressure- less free fall (red middle line). A viable accretion solution with a shock starts at r v and follows any of the black lines until 
it crosses a dotted line, where it jumps to the assumed upstream velocity profile. Above a certain c£f* there is no accreting solution with 
a steady-state shock. The critical value x CTlt , where this happens, is shown with a vertical red dash-dot-dot line for the pressure-less free 
fall upstream of the shock. 



possible, qualitatively similar to the case of Bondi flow 
(blue dotted line) described above. 

Finally, we briefly discuss the "upper" solutions that 
appear when ct is higher than the value correspond- 
ing to the decelerating transonic solution (dashed line 
in Fig. [1]), but still ct < c™ *. We see that the maxi- 
mum shock radius for these solutions is identical to the 
sonic point radius. The shock radius decreases with 
increasing ct until it merges with the normal solution 
branch at c" 1 *. These upp er solutions were found by 
lYamasaki fc Ya mada ( 2005) also in the full problem, but 
these authors show that the upper solutions are unstable 
to perturbations and thus most likely do not occur in a 
real physical system. Because of this, we do not discuss 
the upper solutions further except in Sections 13.21 and 



2.2. The critical sound speed 
In Figure [D we show the critical values rt as a func- 
tion of mass accretion rate M for the three sets of shock 
jump conditions shown in Figure [TJ The three curves 
are similar over a significant range of M despite the rel- 
atively stark differences in the shock jump conditions. 
Specifically, the top and the middle curves closely resem- 



ble the bottom blue line, which is, in fact, the relation 
between the sound speed and the mass-accretion rate 
for is othermal transonic accretion ([Lamers fc Cassinellil 
119991 p. 68). The curvature of this relation in the 
log-log plane of Figure [2] arises from exponential near- 
hydrostatic density profile of the atmosphere, which is a 
good approximation to the decelerating transonic accre- 
tion flow below the sonic point. 

From Figured] it is obvious that for pressure-less free 
fall upstream of the shock the solution described by 14 
does not correspond to any special radial point in the 
flow. However, we also see that the shock position x crlt at 
c™' is constant when expressed in the dimensionless co- 
ordinates since this is the point where the dotted line just 
meets the flow profile. Changes in M, p u or r v only move 
the starting inner point of the flow, but one can always 
find a value of ct = c™*(M, p„, r v ) to get the critical so- 
lution at x clit . By plugging the solution of from the 
shock jump condition in equation (|3Ttj) into equation 

6 We note that this procedure is valid only for determining the 
position of the critical point in isothermal accretion. As seen from 
Figure[T] the critical point occurs whe n the line of all possible shock 
positions (solution to eqs. [2] an d | 3bl . dotted red line) just touches 
the flow profiles (solutions to eq. pQ, solid black lines). Thus, at this 
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Figure 2. The maximum isothermal sound speed cf lt that allows for a shock in the flow at a given M. The PNS parameters — M, rv, 



and color coding — green, red, blue — are the same as in Figure [T] 



we obtain exactly x mt = 3/16 implying that the critical 
condition for existence of steady-state accretion shock in 
isothermal accretion is 



2 = ^= - 1875 ' 



(4) 



where v csc — y/2GM/r is the escape velocity. This result 
is valid when the flow upstream of the shock is in free 
fall and has negligible thermal pressure. We note that 
equation (Q} is valid only at dp*. There are indeed shock 
solutions with x s h oc k > x CIlt that corresp ond to t he up per 
sol ution branch as discussed in Sections 12.11 and 13.21 and 
bv lYamasaki fc Yamadal (|2005[ ). but these solutions still 
have ct < dp*. The actual value of dp* for given M, M, 
r v , and p v is calculated numerically as is demonstrated 
in Figure [2] 

As the ratio in equation (UJ is lower than for the sonic 
point, we call this condition "antesonic" (see footnote 2]) . 
For the case of free fall with non-negligible pressure 
shown with green line in Figure [IJ the numerical fac- 
tor is (5 — v / 21)/2 — 0.2087. Furthermore, we emphasize 
that equation Q does not imply any connection to the 
Parker pointQ or sonic point. Thus, one cannot assume 

point only, the tangents and values are the same, which we then 
use to calculate the position of the critical point, equation JJ). 

7 The Parker point occurs when the right-hand side of equa- 
tion (fTJ vanishes, and in the case of isothermal flow the Parker 
point coincides with the sonic point (Lamcrs & Cassinclli 1999, 
p. 63). 



that the existence of a shock is equivalent to the exis- 
tence of a sonic point i n the context of supernova stalled 
accretion shocks (as in iShigevamall 1 9951 ) . However, if the 
flow upstream of the shock is a Bondi accretion flow with 
the same sound speed, then the critical condition is ex- 
actly the same as the expression for the position of the 
sonic point: (c" 1 * / v esc ) 2 = 0.25. 

Our calculation of the critical condition in equation Q 
also explains why the critical curves for different shock 
jump conditions in Figure [5] are so similar: their values 
can be estimated in the same way as for t he position 
of the sonic point (jLamers fc Cassinellil[l999l p. 68), ex- 
cept that for the shock jump conditions relevant in the 
supernova context the value of x crlt is lower than 0.25 
(eq. 01). 

For the sake of completeness, we repeated the isother- 
mal calculation with pressure-less free fall above the 
shock for a polytropic equation of state P — kp v , where 
k is the normalization factor and F* is the adiabatic index, 
which we vary from 1 to 1.7. We do not assume energy 
conservation through the shock so that k is left as a free 
parameter, which we can vary to obtain its critical val- 
ues. The results are qualitatively similar to those for 
isothermal accretion. For fixed T, there exists a critical 
value of the constant k above which there are no steady- 
state solutions with a bounding shock. As in the isother- 
mal case, this critical value of k corresponds to a critical 
condition that can be written as max(c|/^g SC ) s» 0.19r 
(compare with eq. [4]), where eg is the adiabatic sound 
speed and the "antesonic factor" is linearly proportional 
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Figure 3. Antcsonic factor for polytropic equation of state. For 
each r, we plot max (cr|/t;;? sc ) as a function of k for M ranging 
from -0.01 Mq s -1 (black lines) to —1.0 M® s -1 (red lines). For 
fixed r, the lines overlap. The antesonic factors at each critical 
k are marked with filled circles with color corresponding to M. 
The antesonic factor at critical k is linearly proportional to F 
The horizontal dotted line marks the value for the isothermal case, 
3/16 = 0.1875. 



to the adiabatic index of the flow. We illustrate the de- 
pendence of max (c^/Wg SC ) on T in Figure [3l where we 
plot the antesonic ratio as a function of k for a range 
of M and two values of T. We note that although we 
are able to write down a single local criterion for the an- 
tesonic condition in the adiabatic case, at the radial lo- 
cation where max (c;|/i>g SC ) occurs, this condition should 
be interpreted as a global condition on the flow, and not 
a local one. As in the isothermal case, the physics of 
the critical condition is determined by the ability of the 
flow to satisfy the jump conditions while simultaneously 
satisfying the Euler equations for steady-state accretion. 

The physics of collapse, core bounce and shock stall 
dictate that there is a phase of steady-state accretion flow 
with a standing shock with time-changing M and ct- If 
Ct slowly increases at fixed M, the shock radius advances 
until the critical value c^ 1 * is reached. Physically, c^ 1 ' 
exists because as ct is increased the shock jump con- 
ditions simply cannot be satisfied given the properties 
of the flow upstream and the necessity of maintaining 
conservation of mass and momentum. What happens 
at c" 1 * depends on the properties of the flow upstream 
of the shock. In the case of Bondi flow, the profile can 
continuously transition to a sho ck-less B- t ype s ol ution 
(accre tion breeze ) as shown by iKorevaarl (|1989f ). iVellil 
Ii99il200l and lDel Zanna et all (11998). because there 



is no jump in velocity at c™*. However, in the case rel- 
evant for supernovae, where the flow above the shock 
is in pressure-less free fall, the shock does not degener- 
ate at rt and a catastrophic dynamical transition must 
occur, most likely to a supersonic wind, which can be 
identified with the supernova explosion (iBurrowsl Il987t 
iBurrows fc Goshy| |1993l; iBmrows et~aTlll995fl. A simila r 
conclusion was reached bv lYamasaki fc Yamadal (|2005j ). 
It is likely that time-dependent hydrodynamical insta- 
bilities will modify the exact realization of this scenario. 
Nevertheless, our results from the isothermal accretion 
model are an important step towards understanding the 
mechanism of supernova explosions, and specifically the 
neutrino mechanism, as we detail below. Our isothermal 



model also provides a robust explosion criterion of broad 
applicability (the "antesonic" condition, Section |4~4|) . 

3. STEADY-STATE ACCRETION SHOCK IN 
CORE-COLLAPSE SUPERNOVAE 

While isothermal flow calculations provide insight into 
the physics of steady-state accretion with shocks, we per- 
form more involved calculations to assess the effects of 
individual parameters of the problem for supernova ex- 
plosions. In this Section we first describe our equations 
and boundary conditions and then we proceed by dis- 
cussing how the solution changes when approaching the 
critical neutrino luminosity. We then prove the corre- 
spondence between the critical neutrino luminosity and 
the critical sound speed from Section^ and discuss prop- 
erties of the critical solutions. 

3.1. Numerical setup 
We solve the time-independent Euler equations 



0, 



1 dp 1 dv 2 
p dr v dr r 
dv 1 dP _ 
dr p dr 
de P dp q 
dr p 2 dr v 



GM 



(5) 
(6) 
(7) 



where P is the gas pressure, M is the mass within the 
radius of the neutrinosphcrc r v , e is the internal specific 
energy of the gas, and q = % — C is the net heating 
rate, a difference of heating H and cooling C. These 
hydrodynamic equations couple through the equation of 
state (EOS), P(p,T,Y e ) and e(p,T,Y e ), where T is the 
gas temperature, to the equation for electron fraction Y c 



V — lu e n 
dr 

where l Ve m h+ 



+ L 



„, lj7 cP , and l e - p are reaction rates involv- 
ing neutrons and protons. In Appendix we provide 
equations ([5H7J) in the form of separated radial deriva- 
tives, which is useful for describing the connection be- 
tween the critical luminosity and cSf* of Section [5J We 
also employ, for the first time in the context of calcula- 
tions of the critical luminosity, gray neutrino radiation 
transport using the expression 



— — = -4nr pq, 
dr 



(9) 



where L v is the combined luminosity of electron neutri- 
nos and antineutrinos, which are assumed to be equal ev- 
erywhere. We do not consider neutrinos of other flavors. 
We assume that neutrinos have root-mean-square en- 
ergy of e„_ =13 MeV a nd antineutrinos e Pc = 15.5 MeV 
([Thompson et all 12003). Even this very simple approx- 
imation to radiation transport enables us to assess the 
relative importance of accretion and core neutrino lumi- 
nosity for the supernova explosion. Our choice of equa- 
tion of state, heating and cooli ng function, and rea ction 
rates is a simplified ve rsion of lScheck etaTl J2006) and 
Q 

ian fc Wooslevl (|1996l ). More details are given in Ap- 
pendix [A] 

Five equations ([SH5]) are solved on an interval of radii 
between r v and rs, where r„ is the radius of the neutri- 
nosphere and is fixed, and rs is the radius of the shock, 
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which is left to vary. The matter inside of r v is the PNS 
core and is characterized only by its mass M and lumi- 
nosity L Vt corc and energy of neutrinos and antineutrinos 
that it isotropically emits. 

Wc thus need six boundary conditions to uniquely de- 
termine five functions: p, v, T, Y e and L„, and shock ra- 
dius rs- We first demand that the flow has a fixed mass 
accretion rate M = 4irr 2 pv by applying this constraint at 
the inner boundary. We assume that the neutrino lumi- 
nosity at the inner boundary is equal to the luminosity of 
the PNS core, L v (r v ) = L V COIC . At the outer boundary 
we apply the shock jump condition^ 



pv 



P = P + v 2 s , 

P 1 o 



— V 

2 



■ £ H = -V 



a- 



(10) 
(11) 



where vs — \/Tv CS c is the free fall velocity, and in ac- 
cordance with analy tically estimated accretion data of 
IWooslev et al.l (|2002f ) supernova progenitor models avail- 
able on-lin^3 we choose T = 0.25. The quantity p + is 
the density just upstream of the shock and can be calcu- 
lated from conservation of mass. We also assume that the 
matter entering the shock is composed of iron and hence 
Y c = 26/56 at the outer boundary. The last boundary 
condition comes from requirement that r v is the neutri- 
nosphere for electron neutrinos, which gives a constraint 
on the optical depth 



n Uc p dr = -, 



(12) 



where n Ue is the opacity to electron neutrinos. This con- 
dition is implemented by adding an extra variable corre- 
sponding to the optical depth and requiring that it is zero 
at the inner boundary and 2/3 at the outer boundary. 

The equations and boundary con ditions above are 
solved using a relaxation algorithm ([Press etaDIMl 
p. 753) with logarith mic spacing of grid points (see 
iThompson et all 1200 ll and references therein). The 
number of grid points was usually > 1000. We refer to 
the problem outlined by equations (|5T IT2l) and the physics 
described in Appendix [3] as the "fiducial calculation" . 
As this problem is complex, we occasionally resort to 
various simplifications to make some of our points more 
clear. We describe what exactly we modified with re- 
spect to the fiducial calculation at a various places in the 
text. 

3.2. Properties of the fiducial calculation 

We include neutrino radiation transport and hence our 
solutions are parameterized by the core neutrino lumi- 
nosity L v _ corc and yield neutrino luminosity as a func- 
tion of radius L v (r). We confirm the existence of an up- 
per limit to the core neutrino luminosity £" corc at fixed 



M (jBurrows fc Goshvl [19931: lYamasaki fe Yamadal [20051 
[2006J). We determine L"^ mc by starting the relaxation 
code at a low value of L Vi corc with an initial guess of 
power-law profiles of variables, and then increase L v corc 

8 We do not include nuclear binding energy in our calculations. 
However, the p ertinent effects are qualitati vely assessed in Ap- 
pendix [A] and in Yamasaki & Yamada (2006). 

9 http : //www . stellarevolut ion . org/data . shtml 



in small steps while using the converged solution from 
the previous step as an initial guess for the next step. 
We use bisection to trace L"^ mc to a specified relative 
precision. 

Figure [4] shows profiles of thermodynamic quanti- 
ties for L^, core stepping from to L"^ m . c . Similar 
plots except for L„( r ) were pub lished previously by 
lYamasaki fe Yamadal ([2005L 120061 ). but we reproduce 
them here for the sake of completeness. Our results 
fall within the pa radigm of quasi-steady-state super- 
nova structu re (e.g. iBruennl 1985; M avle fc Wilson|[l988l : 
lJankall200iT ). In particular, we find a steep decrease in 
the density outside r v that flattens when the dominant 
pressure support changes from an ideal gas of nucleons 
to a relativistic gas of electrons and positrons. The tem- 
perature profile is flat at small r, but then drops as the 
relativistic particles become more dominant. Still, the 
temperature changes only from approximately 3 MeV to 
0.5 MeV, much less than the corresponding change by 
several orders of magnitude in density. This implies that 
our isothermal model is a legitimate rough approxima- 
tion of the problem. We discuss in more detail the pro- 
files of the adiabatic sound speed c?, total specific energy 
in the form of Bernoulli integraFI 



P GM 

£+ , 

p r 



(13) 



and neutrino luminosity L v that are to our knowledge 
not previously discussed in the context of the calculation 

rcrit 
01 -^corc- 

The total specific energy 3$ is an increasing function 
of L„, corc in the whole region of interest. The maximum 
of SS occurs at r ga i n , which is defined as g(r ga i n ) = 0. 
Both r ga j n and ^(r ga ; n ) increase with increasing L v cmc . 
We note that even at £™c 0re the energy is significantly 
negative in the entire region of interest. It is certainly 
not the case that positive 5$ correspo nds to explosion 
(L ViCOle = Lf" c t D1 . c : lBu7ro~ws et al.l[l99m . 

The ad iaba tic sound speed eg was calculated using 
equation (|A9|) and is plotted in the bottom left of Fig- 
ure 2] as a ratio to the local escape velocity v csc (r) . Sim- 
ilar to we see that an increase of L„, corc results in a 
net increase of c$ between r„ and rs- The maximum of 
cs occurs at gradually larger radii for increasing L UyCorc , 
similar to the gain radius r ga ; n . However, the maximum 
of c|/Ug SC does not exactly coincide with r ga j n ; in fact, 
'"gain is always about 1 to 2% larger for the particular 
calculation shown in Figure [4j 

The middle center panel of Figure [4] shows profiles of 
L v (r). As expected, L v grows in the cooling layer as a 
result of net neutrino and antineutrino production and 
decreases in the gain layer as a result of heating. How 
much does the accreting material cooling via neutrino 
emission contribute to the total neutrino flux emanating 
from below the shock? We define the accretion neutrino 
luminosity Lf, cc as the difference between the neutrino 
luminosity at the shock and the luminosity of the core, 

Ll cc = L„(r s )-L^ col . c . (14) 

The relative contribution of the accretion luminosity 

10 The importance of individual terms contributing to 8$ is as- 
sessed in Appendix |B1 
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Figure 4. Variables of interest as the core neutrino luminosity L„, CO r C approaches the critical value i"c 0re for M - 
M = 1.4 Mq and r v = 50km. Profiles of density p, velocity v, temperature T, electron fraction Y c , neutrino luminosity L h 



-0.5MQS" 1 , 
net heating q, 



adiabatic sound speed eg, internal specific energy e, and total specific energy , 
The profiles at L^ T1 ^ OTC are shown with blue thick lines. 



' for < L v 



: 6.41 X 10 52 ergs s x are shown. 



to the total neutrino luminosity decreases with grow- 
ing L^coro- For core = 0, the accretion luminosity 
is the sole contributor to the total luminosity of about 
1.4 x 10 52 ergs s _1 , but for the critical core luminosity of 
6.41 x 10 52 ergs s _1 the accretion luminosity contributes 
28% to the total neutrino output. Clearly, the accretion 
luminosity is a sub-dominant, but important part of the 
total neutrino output. In Section [3.41 we investigate how 
including the neutrino radiation transport and hence ac- 
cretion luminosity affects L™^ for a range of M. 

We now turn our attention to another aspect of steady- 
state accretion shocks: the existence of two solutions 
with differ ent shock radii for some valu es of L v ^ core , dis- 
covered bv lYamasaki &: Yamadal (|2005h . These authors 
also performed linear stability analysis and found that 
the upper solution branch with higher rs is unstable, 
the lower branch is stable, and both branches meet at a 
marginally stable point at L™* orc . We confirm the pres- 
ence of the upper solution branch with our code. We 
switch to the upper solution branch by first calculating 
a sequence of models on the lower branch up to L"^ mc . 
We take a solution close to L£"* ore , manually scale up the 
radial coordinates by about 50% and let the relaxation 
algorithm converge to a solution on the upper branch. 

We find that the profiles of thermodynamic variables 
for the upper solution branch are not fundamentally dif- 
ferent from those of the lower branch presented in Fig- 



ure HI and hence we do not present them here. We plot 
in Figure [5] the shock radii and total specific energies 
at several radii of interest for both branches, and as a 
function of L v , core for several values of M . We see that 
the upper branch of solutions extends to larger shock 
radii, but ends at a finite radius and non-zero L Ut coro 
(left panel, Figure [5]). We specifically checked that this 
is not an artifact in our calculation. We note however, 
that for different input physics the upper branch can 
end at much larger radii, consistent with infinity in our 
code. We also find that the specific energy of the upper 
solution is always higher than that of the lower solu- 
tion for the same L v ^ coro at all radii. In analogy with 
other physical systems that have two states with differ- 
ent energy, we can argue that the higher energy upper 
soluti on branch is unstabl e to tr ansition, in agreement 
with lYamasak i fc Yamadal (|2005|) who find that the up- 
per branch is indeed unstable to radial perturbations. 
We do not discuss the properties of the upper branch in 
greater detail, because it is unlikely that these solutions 
would be realized in a real physical system. However, 
the existence of the upper solution branch is analoguous 
to the upper solutions in Figure [TJ 

3.3. Correspondence with isothermal accretion 

It has been clearly established in ID calculations that 
there exists an upper limit to the core neutrino hum- 
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Figure 5. Properties of the two branches of solutions for M 
(green) . Left: Shock radius rg as a function of Lv, < 
mass accretion rates. The black dotted line denotes r v . The solid and dashed parts of each line denote the lower (stable) and upper 
(unstable) solution branches, respectively. Right: Total specific energy 33 at r„ (bottom lines), r ga i n (top lines) and rg (middle lines). 
Color coding is the same as in the left panel. 



nosity L"^ OIC that allows for a steady-state accretion 
shock (fBurrows k. Goshvl Il993t lYamasaki fc Yamadal 
2005, 2006), and this limit has been confirmed to exist 
in multi-dimensional dynamical supernoya sim ulations 
(|Murphv fc Burrowd[200llNordhaus et alJl2010D, and in 
simulations of the stalled a ccretion shock ( Ohnish i~et al.l 
I2006t llwakami etldl 120081 ). To illustrate how the up- 
per limit on the core neutrino luminosity is related to 
the limit on the sound speed in isothermal accretion dis- 
cussed in Section [5J we consider the total specific energy 
of the system 38, which has the property that its radial 
derivative is proportional to the net deposited heating 
and cooling: 

dSS 4irr 2 p . , . 

q. (15) 



dr 



M 



It is well established (e.g. Qanka 2001) and we confirm 
this in Appendix [B] that the contribution of the kinetic 
term v 2 /2 to the total specific energy budget in equa- 
tion (fT5|) is negligible. At a fixed radius r in the gain layer 
where q w L u ^ COIC /r 2 , the increase of L Vi cor c will result 
in an increase of enthalpy h — e + P/p, and consequently, 
of the adiabatic sound speed, as h ~ c s . We showed in 
Section HOI that cs as well as 3$ at a given radius are in- 
creasing functions of L Uy core everywhere, not only in the 
gain layer. This correspondence between L„, corc and cs 
suggests that L£ rl * ore is a realization of c™* (Section [5]) 
in the context of the neutrino heating mechanism. 

As shown in Section [51 c" 1 * occurs when the flow can- 
not satisfy the shock jump conditions at any point i n 
the accretion flow (see also lYamasaki fc Yamadal 12005). 
We now prove that the same mechanism is responsi- 
ble for the existence of L"^ oie in the supernova con- 
text. Essentially, we are trying to construct an equiv- 
alent of Figure [T] for non-isothermal flows and more real- 
istic physics. In order to make this feasible we consider a 
simplification of the problem by setting dY c /dr = and 
dL v /dr = (eqs. [8] and [9]) and setting the heating 
and cooling terms to % = 1.2 x 10 _18 cm 2 g~ 1 L V:Cmc /r 2 
and C = 2 x 10 18 (T/MeV) 6 ergs s _1 g"\ respectively. 
We also modify the boundary conditions by requiring 
p{t v ) = 3 x 10 10 g cm -3 instead of fixing r v . With 



these changes, the problem simplifies considerably, be- 
cause p(r v ) and v(rv) are now fixed and the only vari- 
ables left are T(r v ) and r$, which are determined by the 
two shock jump conditions. Consequently, if we know 
the value of T(r y ), the flow profile is now uniquely de- 
termined as an initial value problem. 

We want to show how the solution for the flow struc- 
ture with a shock relates to all possible accretion solu- 
tions. For a given L v _ corc and T(r u ) we calculate the 
A- type solutions (those reaching = — 1 in Figure [T]) 
by setting v 2 — c 2 s at the outer boundary, which gives 
a continuum of solutions for a range of T{r y ). Solutions 
going through the sonic point are obtained by setting 
the numerator and denominator of the momentum equa- 
tion (IA2[) simultaneously to zero at the outer boundary. 
The accretion breezes (B-type solutions in Figure [T]) are 
obtained from the A-type solutions by setting v to a spe- 
cific fraction of c$ at a fixed radius corresponding to 
v = cs in an A-type solution. By varying this fraction, 
we get a continuum of flows with different T(r„). We 
extend the breezes to larger radii with an initial-value 
integrator of stiff equations based on algorithms from 
I Press et all (|1992h and using the results from our relax- 
ation algorithm as a starting point. We verified that for 
B-type solutions (subsonic accretion breezes) the ratio 
|i>/cs| is always less than unity. 

In the left panel of Figure|B]we show the position of the 
shock radius rs (solid) and the sonic point r son i c (dashed) 
as a function L v ^ coro . In the right panel, we show the val- 
ues of T(r u ) for the same solutions as in the left panel. 
We see that r son j c exhibits qualitatively similar behavior 
to rg , except that for low values of L u , coro there is no 
sonic point with a finite r son j C . Furthermore, the posi- 
tion of the shock is well separated from the sonic point, 
except for the tip of the upper solution branch, where 
these two appear to merge, similar to what happens in 
isothermal accretion in Figure [TJ The same trend is seen 
in T(r v ) (the only free parameter of the problem for fixed 
L u , core), which suggests that the upper branches of shock 
and sonic point indeed coincide at some L„ iCorc . 

In order to explore the situation in more detail, we 
selected four values of L Ut core marked by vertical dot- 
dashed lines to illustrate the four configurations visible 



10 



Pejcha & Thompson 




- 1(1) 

- j 


(2) 

B 


(3) (4) 


- 


- B y 


/ 

/ 






- j / / 
i ••' / / 
~A A / ' 

' \ / 1 
i / / 
i / / 


B 






- i/ ' 

! A / 

/ 

! 




1 





4.0 



3.6 % 



3.4^ 



3.2 
3.0 



50 



100 



ergs 



-l 



Figure 6. Overview of the flow structure in the simplified supernova setup with L'^J 1 ^ OTe = 9.07 X 10 52 ergs s 1 . Vertical dash-dotted lines 
denoted by numbers mark values of L^, core selected for detailed study in Figure [7] Left: Shock radius rg (red solid line) and sonic radius 
r sonic (black dashed line) as a function of L VyCorc . Right: Temperature at the neutrinosphcrc T(r v ) for the solutions with a shock (red 
solid line) and with a sonic point (black dashed line). Regions of different solution types are marked with letters A and B, as in Figure [T1 
The dotted line is the separatrix between A- and B-type solutions where the sonic point lies at infinity. 



in Figure [BJ (1) low L VjCOTe with one solution with finite 
shock radius and no finite sonic radius, (2) intermediate 
L v core where two solutions are present for both the shock 
and the sonic point, (3) L Vi corc > L"^ ore but with two 
sonic point solutions still present, and (4) high L v corc 
with no solution for either the shock or the sonic point. 
In Figure [7] we plot flow solutions for each of these lu- 
minosities (black solid lines) with shock jump conditions 
overplotted (green solid lines). The Figure should be in- 
terpreted in a way similar to Figure [U We found that 
plotting the Mach number of the flows only obfuscates 
the issue because the lines cross in a complicated manner; 
instead we plot momentum pv 2 + P and specific energy 
38. In this parameter space the different flow lines gener- 
ally do not intersect and the shock jump conditions are 
a function of only the radius as can be seen from the 
equations (fTU)) and (|11[) . 

The upper left part of Figure [7J shows flow profiles for 
low core , where there is only one solution for the shock 
radius (far left vertical dash-dotted line in Figure [6]). 
Each solid line corresponds to a different value of T(r v ). 
For low T(r„), we find that the flows are of type A. That 
is, they reach a point where v 2 = c| (not a sonic point; 
shown with black dots in Figure [Jj. For increasing T(r„), 
the intersection of the flow with the shock jump condi- 
tions moves to larger radii. A-type solutions extend to 
T(r v ) ~ 3.32 MeV, where the point where v 2 — c| moves 
to infinite radius. For yet larger T(r„) we are able to find 
only B-type solutions — subsonic accretion breezes - 
with |v/cs| < 1 everywhere. For a single value of T(r v ) 
the profiles of momentum and energy of the flow inter- 
sect with the shock jump conditions at the same radius, 
which means that only at this T(r v ) is a solution with a 
shock possible (red line) . The vertical dotted lines mark 
radii of intersection of momentum of the flow with the 
momentum shock jump condition. We then plot these 
radii in the energy panel and we see that for flows lying 



below the red line the intersection of 38 with the ap- 
propriate shock jump condition occurs at systematically 
smaller radii than what is required by the momentum 
profiles. Conversely, for curves above the red line the 
intersection occurs at larger radii. Thus, the momentum 
and energy at the shock can only be conserved simulta- 
neously for a single T(r u ) at L VlCOIe = 5x 10 51 ergs s~ x , 
which then defines not only T(r !/ ), but also rs- 

The upper right panel of Figure [7J displays flow profiles 
for intermediate L u ^ COTC , which allows for two solutions 
for the shock and the sonic point (vertical dash-dotted 
line labelled (2) in Figure [BJ). The topology of the solu- 
tions for this -L„ iC orc differs from the previous case. The 
A-type solutions exist only for a small range of temper- 
atures (3.517MeV < T(r v ) < 3.675 MeV). There are 
two transonic solutions going through two distinct sonic 
points at finite but vastly different radii, which separate 
A-type solutions from the B-type solutions, (which occur 
for T{r v ) < 3.517MeV and T{r v ) > 3.675 MeV). Due to 
severe stiffness of the equations, we were not able to fol- 
low the B-type solutions to very large radii. In contrast 
to the upper left panel of Figure [7J at lower L Vi corc , the 
flow profiles intersect the shock jump conditions at just 
two inner boundary temperatures shown by the two solid 
red lines. The upper solutions are analoguous to those 
obtained in the isothermal model in Figure [T] 

The lower left part of Figure [7] explores the region 
with -L„ iC ore > £"corc where two solutions for the sonic 
point are still present. The topology of the solutions is 
similar to the previously discussed case with L Ut corc = 
5 x 10 52 ergs s _1 , except that the A-type solutions exist 
for an even smaller range of inner boundary temperatures 
(3.852 MeV < T(r u ) < 3.875 MeV). However, in this 
case, there are no simultaneous intersections of the mo- 
mentum and energy profiles of the flows with the shock 
jump conditions. This is explicitly shown by the verti- 
cal dotted lines, which mark several radii of intersection 
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Figure 7. Detailed study of the structure of flows in the simplified supernova setup summarized in Figure [6] The four panels show 
flow structure at different values of L v , cove corresponding to different situations: single shock radius and no sonic point (upper left), two 
shock radii and two sonic-point solutions (upper right), no shock solution and two sonic-point solutions (lower left), and neither shock or 
sonic-point solution (lower right), fn each panel the flows that satisfy the shock jump conditions are shown with red solid lines, flows going 
through the sonic point (filled triangle) are shown with black dashed lines, and all remaining flows are shown with solid black lines. A-typc 
solutions going through v 2 = cj| end with filled circles. Regions of different solution types are marked with letters A and B, as in Figure l6l 
Green solid lines mark the shock jump conditions. Additionally, the dotted vertical lines in the upper left and lower left panels mark radii 
of intersection of the flows with the momentum shock jump condition. 
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of the momentum profiles of the flows with the momen- 
tum shock jump condition. Clearly, the intersections of 
the energy profiles with the energy shock jump condi- 
tions occur at systematically smaller radii. This means 
that at this L v ^ coro no flow profile can connect to the 
free-falling matter upstream, because it is not possible 
to simultaneously satisfy conservation of mass, momen- 
tum and energy at the shock. Although the mismatch 
in radii decreases as we go to the higher T{r u ) of B-type 
solutions, they eventually cease to cross the energy shock 
jump condition altogether. Overall, as implied by Fig- 
ure El for L^corc = 9.5 x 10 52 ergs s _1 it is impossible 
to find a flow solution that simultaneously conserves mo- 
mentum and energy at the shock. 

For the sake of completeness we show in the lower right 
panel of Figure [7] the flow profiles for high L Ut corc , where 
there is no possible solution with a shock or a sonic point. 
All flows displayed are B-type solutions with |u/cs| < 1. 
Despite the fact that we were not able to extend some of 
the solutions to large radii due to stiffness of the equa- 
tions, v/cs of the solutions systematically converged to 
zero at large radii. 

With the L UtCme -T{r v ) space systematically explored 
we can now label regions of A and B solutions in Fig- 
ure El in analogy with Figure Q] We note that the more 
physically realistic simulation displayed in Figure [5] has 
exactly the same behavior except that sonic points ex- 
ist at finite radii all the way to L v _ corc = 0. The upper 
shock solution branch, however, ends as indicated in Fig- 
ure [5j As can be seen from Figure El the position of the 
shock when L Ut corc = L"^ me does not align with the po- 
sition of the sonic point and hence a smooth transition 
to a B-type solu t ion is not possible , unl ike the results of 
IKorevaarl (pSflh . rVeliil (fl9M 12001 and lDel Zanna et al.l 
(1998))! With the same logic as in Section [2] we ar- 
gue that at ££ rl c 0re the structure of the accretion flow 
catastrophically changes: if L U COIC > L£ rl * ore the shock 
moves outwards dynamically and a neutrino-driven wind 
with a sonic point below the shock is established. How- 
ever, multi-dimensional time-dependent hydrodynamical 
instabilities might modify the realization of this process. 

3.4. Properties of the critical solutions 

In Figure we present the critical core neutrino lumi- 
nosity L™^ te versus the mass accretion rate, the critical 
curve, as calculated with the fiducial physical model as 
described in Section 13.11 Before we dive into properties 
of the critical curve, we investigate whether it is a robust 
phenomenon from the point of view of input physics. To 
this end, we simplify the problem by setting the electron 
chemical potential to zero, shutting off neutrino radia- 
tion transport (dL^/dr = 0), electron fraction evolution 
(dY c /dr = 0), and choosing simplified heating and cool- 
ing prescription with H s = 1.2 x 10~ 18 L u ^ COIC /r 2 and 
C s = 8.19 x 1(T 43 T 6 , respectively. The critical curve 
for this modification is shown in Figure [8] with a dashed 
black line. With respect to the simplified version of the 
calculation we make additional individual changes to see 
how the critical curve changes: we modify cooling by 
multiplying it by a factor of 5, making it proportional to 
T 3 instead of T 6 , or shutting it off altogether. We mod- 
ified heating by making it proportional to r -3 , constant 
at all radii or adding to H s heating motivated by Alfven 



Thompson 

wave energy deposition (|Metzger et al.l [2007). We also 
changed the EOS to be that of a pure ideal gas. We also 
tried an EOS of pure relativistic particles, but we were 
able to converge to solution only by relaxing the condi- 
tion on optical depth and replacing it with fixed density 
at r v . We thus do not show the critical curve here. 

Figure [5] shows that even drastic changes to the func- 
tional form of the heating and cooling functions and EOS 
consistently yield critical values of L v ^ coro . Naturally, the 
exact value of L° n ^ OIe at a given M changes considerably. 
However, it is interesting that we do not see wildly dis- 
crepant slopes of individual lines in Figure [8] Specifi- 
cally, changing the form of the heating affects the slope 
only very slightly, while the form of cooling has a much 
more pronounced effect, except when simply multiplied 
by a constant factor. This can be explained by the fact 
that cooling couples directly to the temperature, and the 
temperature gradient is directly responsible for the hy- 
drostatic equilibrium of the region when relativistic gas 
dominates. 

In order to obtain dependencies of i"c re on * ne P a ~ 
rameters for the fiducial problem (solid line in Figure H]) , 
we determine critical core luminosities on a grid of pa- 
rameter values chosen as M g {1.2, 1.4, 1.6, 1.8, 2.0} M©, 
r v e {20, 40, 60} km, t v e {1/6, 1/3,2/3}, and 30 values 
of M logarithmically spaced between —0.01 Mq s _1 and 
—2 Mq s _1 . Although in principle the neutrinosphere 
radius r v and core luminosity L Ul corc are connected by 
physics of neutrino diffusion out of the PNS, we consider 
them separately here. The critical curves for sample val- 
ues of parameters are shown in the left panel of Figure [9j 
In addition to the well-known dependence of L"^ OT0 on 

M, we see that L™* ore increases with increasing M and 
decreasing r„ and t v . Lower values of t v also increase the 
curvature of the lines in the M-L"' 1 * ore space. In order 
to quantify these observations we fit to all our results a 
combination of power laws to get 

£^ = 8.18 x 10 53 ergs ^ {^j ^ * 

The relative differences between calculated and fitted 
values of L c J^ mc are plotted in the right panel of Fig- 
ure [9l On the first sight, the multiple power-law fit of 
equation (|16j) is not perfect. However, for most of the 
data the power- law fit is good within 30% of L^ OTe . 
There are systematic trends in the residuals as a func- 
tion of M in the sense that all lines are curved upward: 
at low and high M the value of L^ orc is slightly higher 
than what would correspond to the mean power law fit- 
ted to the data. We can understand this on the basis 
of our isothermal accretion model and its critical curve 
presented in Figure[2j which also exhibits an upward cur- 
vature of the critical curves, which is a manifestation of 
the nearly-hydrostatic exponential density structure of 
the flow (see Section [2]). 

In the fit of equation (1161) we did not include the depen- 
dence on the remaining parameters of our calculation, in- 
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Figure 8. Effects of different physics on the critical curve for M = 1.4Mq, r„ = 50km. The solid black line shows the fiducial model, 
while the black dashed line denotes the model with simplified physics (see text for details), which was then subject to additional changes. 
We consider no cooling, C = (solid blue), cooling increased by a factor of 5, C = 5C S (blue dashed), cooling proportional to T 3 , C = 
8.19 X 10 -11 T 3 (blue dotted), heating proportional to r -3 , H = H s X (r„/r) (red dashed), constant heating rate T-L = 1.2 X 10~ 18 L„, corc/r 2 . 
(solid red), neutrino and Alfven waves heating, W = "Hs + 10 21 ergs s _1 g -1 (r^/r) 2 exp[(l — r/r^)/5] (dash-dotted red; Mctzgcr ct al. 2007), 
and EOS in the form of pure ideal gas (solid green). 




Figure 9. Left: Curves of critical core neutrino luminosity as a function of M. The black solid line is for M = 1.6 Mq, r„ = 40 km and 
T u = 2/3. Red dotted lines are for the same parameters except M = 2.0 Mq (upper line) or 1.2 Mq (lower line). Green dashed lines are 
for r„ = 20 km (upper line) or 60km (lower line). Blue dash-dotted lines are for t„ = 1/3 (lower line) or 1/6 (upper line). Right: Relative 
residuals from the fitting of equation d 1 6 1) to the ensemble of calculated critical luminosities. Each line corresponds to a sequence of M 
from —0.01 to —2Mq s _1 with other parameters fixed. 
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Figure 10. Critical shock radii relative to r u as a function of 
^if core f° r calculations presented in Figure [9] Different colors cor- 
respond to different values of M . 

eluding the incoming electron fraction Y c (rs), and energy 
of the neutrinos e„ c and antineutrinos ep c . We performed 
a simple check to see how changes in these parameters 
affect L£ rl ' orG . We find that changing Y e (rs) from 0.45 
to 0.5 shifts the critical luminosity by at most 1%. As 
suggested by the physics of the neutrino heating, £"'c 0rc 
scales as e~ 2 . We confirm this scaling to within 2% by 
varying e Vc from 8MeV to 14MeV while keeping the ra- 
tio e^/e^ fixed at 13/15.5 ~ 0.839. 

One of the uncertainties in core-collapse supernova 
simulations and neutron star physics in general is the 
EOS of dense matter. For example, a softer high- 
density nuclear EOS generically leads to smaller PNS 
radii at earlier times after collapse and bounce, but also 
leads to higher average neutrino energies, which mig ht 
be favorable for expl o sion (e.g. iBaron et al.l Il985allbl; 
Sumivoshi et al.l 120051 : Uanka et al.l l2005t iMarek et al I 
20091 ). Can we put some constraints on the EOS us- 
ing our results on L"^ OIC l In more realistic simulations, 
L v , corej r v and the temperature of emitted neutrinos 
T v ~ e„ e are connected by the physics of diffusion. In 
the simplest case, we can assume that these quantities 
follow the black-body law 



4L, 



7ircrT,l ' 



(17) 



where a is the Stefan-Boltzmann constant. In addition 
to equation (TTcj)) we found that ££ rl * ore ~ e^ 2 to a very 
good precision. Eliminating the dependence on e Ve using 
equation (1171) . we find 



LI 



oc t,7°- 14 M 



1.23 



(18) 



Therefore, larger neutrinosphere radii suggest lower 
^ifcore: when all other quantities are equal. Thus, stiff 
equations of state that yield higher r„ are potentially 
preferred. We note that this result is enabled by su- 
perlinear dependence of L c J^ mc on r v in equation (|T6)) . 
On the other hand, soft EOS leads to a faster con- 
traction of the PNS, higher core luminosities and neu- 
trino en ergies, which can be favorable for the explo- 
sion too (IScheck et al.ll2006l [20081: IMarek fc Jankall200a 
IMarek et al.N2009h . 

Finally, we briefly mention the behavior of shock radii 



as r§ rit /Vr, oc (L£ 

T crit 
v. core 



at L 



crit 



§ rlt . In Figure [TUl we plot r§ rlt jr v as a func 



tion of i^core f° r data from Figure [HI We see that 
the points lie close to a line described approximately 
)-0.26 rpjjg wea k dependence on 

can probably be attributed to changes in the ef- 
ficiency of cooling as a function L£, rl * ore an( i hence M 
(see Section 15.31 for a related discussion) . There are no 
noticeable residuals as a function of r v , M, and t v . 

3.5. Interpretation of L" l ^ mc and the importance of L^ cc 

As we have demonstrated in Section 13.41 the critical 
luminosity depends on M, M and r v . A natural combi- 
nation of these variables to have the dimension of lumi- 
nosity is GMM/r u> an expression for power released by 
accretion of material. Although the powers of individual 
quantities are somewhat different than those in equa- 
tion (|16p. one might wonder whether these expressions 
for the critical luminosity are intimately related in the 
sense that -£™c 0r e i s directly proportional to L^ cc . This 
interpretation is especially tempting, because dynamical 
supernova models seem to be on the verge of explosion 
over a wide range of M, which would suggest that a sig- 
nificant fraction of the neutrino luminosity is supplied by 
cooling of the accreting gas. Any change in M yields a 
proportional change of Lff c , which again somehow pro- 
portionally changes the critical luminosity while keeping 
it higher than Lf, cc . This picture is, however, incorrect. 

To illustrate why, we plot in Figure [TTJ L"^ OTC along 
with the neutrino luminosity leaving the shock L v (rs) 
at L"^ m . c , and their difference, L^ cc . We also plot 

GMM/r u , corresponding to the maximum achievable 
^acc £ or a CQ ^ accr etion flow. We see that L^ cc is al- 
ways a small fraction of -L£ rl * ore , and that even the max- 
imum possible value of Lf, cc falls short of the neces- 
sary luminosity by a fac tor of several as found also by 
iBurrows fc Goshvl (|l993n . Clearly, the major source of 
neutrinos powering the explosion is the core of the PNS. 
However, the accretion luminosity is not entirely unim- 
portant. We plot in Figure [TT] L"% Ie for a calculation 
with the same parameters but with neutrino radiation 
transport shut off. Including the radiation transport in 
the calculation has a positive effect of lowering the criti- 
cal core luminosity by about 8% at M = —0.01 M§ s _1 
and by about 23% at M = — 2M Q s _1 , respectively. 

Instead of regarding equation (|16[) as a manifestation 
of accretion luminosity, we can understand it in the terms 
of the critical condition on the sound speed in isother- 
mal accretion given in equation @. The dependencies 
of equation (|16|) can be then explained in simple terms. 
Increasing the mass of the PNS core makes the escape 
velocity higher and velocity profile steeper, and more 
heating is then necessary to get a mild enough velocity 
profile that it can no longer connect to the shock jump 
conditions (Section [5]). Increasing r„ shifts the region 
of interest to weaker gravitational potential with lower 
escape velocity and hence lower neutrino luminosity is 
required to get a mild velocity profile that again does 
not connect to the shock jump conditions. Decreasing t v 
with all other parameters fixed means less neutrinos ab- 
sorbed and clearly a higher critical neutrino luminosity 
is required. 

To illustrate this in greater detail, we can write 
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Figure 11. Effect of radiation transport and accretion luminos- 
ity on the critical curve. The black solid line shows critical core 
neutrino luminosity L^, T1 ^ OTe for M = 1.6 Mq and r„ = 40 km, the 
black dashed line is the neutrino luminosity at the shock L„(rg), 
and the black dotted line is their difference, the accretion luminos- 
ity L J cc . The blue dash-and-dotted line is the maximum accretion 
luminosity, GM M /r„ . The red solid line is a critical curve for the 
same parameters except that neutrino radiation transport was shut 
off. 

q/v ~ deg/dr ~ c|/rv and hence c| ~ h ~ qr u /v ~ 
p v r v L v> core/ where /i is the enthalpy and p,, is the den- 
sity at the neutrinosphere. Plugging this expression as 



dp* in equation (j4}, we obtain L™ 



crit 
core 



GMM/{ Pv rl). 



The product p v r v is proportional to the optical depth 
t„, because most of opacity to the neutrinos is concen- 
trated close to the neutrinosphere. Thus, by relating c" 1 * 
of the isothermal accretion model through Euler equa- 



tions to L c . 



we arrive through a simple estimate to 



essentially the same expression for L c J^ mc as in the di- 
mensional analysis in the beginning of this Section. The 
power-law dependence of these estimates is of course dif- 
ferent than in the empirical fit in equation (|16p due to 
nonlinear nature of the problem. 

Several authors have produced critical curves based 
on either steady-state calculations like our own, or hy- 
drody namic simulations. Thro ugh steady-state calcula- 
tions, iBurrows fc Goshvl (|1993f ) obtained a critical curve 
that is best fit by a power law, L C J% TC ~ M 0A3 . This 
power-law slope is lower by about 0.3 than what we 
mea sure. This differen c e can be reconciled by realizing 
that IBurrows fc Goshvl ([19931 ) assumed that the neutri- 
nosphere radius r„ is determined by L Vt corc and neutrino 
temperature T v through a black-body law. If we substi- 
tute equation (|17|) to equation (1161) to eliminate r u , we 
get 

(19) 



- crit 



-0.114 7i f l. 02 ji V0. 401,0. 68 



The power- law slope in M is about 0.40, quite close to 
IBurrows fc Goshvl (fl99"l . 

In Fig ure [T2"] we compare ou r results to the critical 
curves of iNordhaus et all (|2010f) based on hydrodynamic 
simulations. In order to make our calculations compat- 
ible with theirs, we couple r v to L v corc through equa- 
tion (HU) with T v = 4.5 MeV. We keep M fixed at 1.4 M , 
which makes the s lope of our crit i cal cu rves slightly 
milder than that of INordhaus et all (|2010f ) : for a given 
progenitor model lower, but fixed, L v _ corc means that the 
explosion occurred later when M is higher, because the 



PNS had a longer time to gain mass, and hence L"^ OIC 
is hi gher. We see tha t our r esults are in good agreement 
with INordhaus et al.l (|2010[) and that including heating 
from accretion luminosity reduces L c J l ^ orc by an amount 
comparable to going from ID to 2D or from 2D to 3D. 

4. EXPLOSION CONDITIONS 

Our goal in this section is to find out how well some 
of the conditions proposed in the supernova literature 
for reviving the stalled accretion shock work in diagnos- 
ing £™* ore within the context of steady-state calcula- 



tion. Similar work was performed bv lMurphv fc Burrows! 
(2008) with dynamical simulations. Our code does not 
allow for time evolution, but we can directly and exactly 
compare how a given explosion condition relates to the 
derived steady-state value of £™core at various values of 
M. 



4.1. Positive acceleration condit ion of \Bethe & Wilsorl 

(MM) 

IBethe fc Wils^nl (Il985h were the first to analyze the 
delayed neutrino mechanism. They argued that shock 
recession creates a steeper pressure profile that in turn 
overcomes gravitation and gives an outward accelera- 
tion, r > 0. In the language of our steady-state model, 
r = v(dv/dr) > is the statement of their critical ex- 
plosion condition. From Figure [4] we clearly see that 
below the shock the flow decelerates everywhere, except 
maybe for a small range of radii, where v is basically 
constant. We also find through an investigation of the 
velocity profiles that having dv/dr = somewhere in 
the flow does not correlate with approaching the criti- 
cal curve at different M. In the discussion of our toy 
model in Section [5] and in Appendix [B] we show that 
velocity and acceleration terms are subdominant in the 
momentum and energy equation of the problem and that 
£"core does not depend on these terms. From this ev- 
idence w e conclude that t he va nishing acceleration cri- 
terion of IBethe fc Wilson] (|1985l ) does not coincide with 
the critical neutrino luminosity. 



4.2. Explosion condition of \Jankd i200A) 

Uankal ()2001f ) presented a quasi-time-dependent ana- 
lytic analysis of the conditions necessary for revival of a 
stalled accretion shock. He found that the shock expe- 
riences both expansion and outward acceleration when 
both the mass and energy in the gain region grow with 
time. These conditions define two critical lines in the 
L u , core~M plane that enclos e a reg i on fav orable for ex- 
plosion. As the appr oach of Uankal (120011) is c ompletely 
different from that of IBurrows fc Goshvl ( 19931 ) and also 
ours, his critical condition is also very different . As 
can be seen in Figures 4 and 5 of ITfankal (pOOlh . the 
critical luminosities are decreasing function of increas- 
ing |M|, unlike the plots in IBurrows fc Goshvl (fl993l) . 
lYamasaki fc Yamadal IpOOH [2006), and our Figure [Til 
where the critical luminosities increase with increasing 
\M\. 

In Figure we plot as a function of L u>COTe and 
M the mass in the gain layer, M Ka ; n = f rs Anr 2 pdr, 
total mass between the neutrinosphere and the shock, 
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Figure 12. Comparison of our calculations with the results of Nordhaus ctal. ( 2010) (red, green and blue lines with points) . The black line 
is our fiducial model with M = 1.4 Mq and r v calculated from equation H17I I, while the dashed line is for our calculation with dL„/dr = 

by an amount comparable to going from ID to 2D or from 2D to 3D. 



(eq. [9]). Heating from the accretion luminosity reduces L" 

Mtot, energy in the gain layer, £; gain = JJ^ Airr 2 p(e - 

GM/r)dr, and the total energy between the neutri- 
nosphere and the shock E to t- We see that for a model 
with fixed mass accretion rate the masses grow with in- 
creasing L„, coro as £"core i s approached. While the total 
energy at any specific radius in the gain layer increases 
with growing L v> corc (see Figure @] for a plot of the re- 
lated Bernoulli integral), the total energy within that 
region decreases, because the gain layer gets physically 

larger with increasing L v COI< >. 

The critical conditions of Janka (2001) demand that 
the mass and energy in the gain region increase with 
time. As there is no time dependence in our steady-state 
models and thus this condition is not directly testable, 
we might emulate the time dependence with a sequence 
of steady-state models. If we fix M and increase core, 
following either of the lines in Figure Q21 M ga i n grows 
and -Egain decreases. On the other hand, if we fix L u , corc 

to, say, 5 x 10 51 ergs s _1 and decrease \M\, which cor- 
responds to moving vertically from red through blue to 
black lines in Figure 1131 M ga i n first drops and then in- 
creases. Similarly, E ga \ n first increases and then just be- 
fore reaching the critical curve it decreases. Thus, from 
the poi nt of vi e w of o ur steady-state models, the condi- 
tions of Danka (200j]) are not equivalent to the critical 
neutrino luminosity. 

4.3. Advection time vs. heating time 



The ratio of advection to heating time £ a dvAheat has 
been extens ively used as a metric of hydrodynamic simu- 
lations (e.g. IThomp son 2000: IThompson & Murravll2001l: 
Thompson et al.l 120051: iBuras et al.l l2006bt iScheck et al.l 
20081) . iMurphv fc Burrows! (|2008t) found that this condi- 



tion is only a rough diagnostic of the critical condition 
for explosion. 

Here, we specifically test two definitions of advection 
and heating times. The first is an integral cond ition on 
the gain layer due to IMurphv fc Burrows! (|2008D 



adv 



dr 



heat 



[ rs 47rr 2 pedr 
J^ S Airr 2 pqdr 



(20a) 



(21a) 



We plot in Figure Q3] values of the ratio t 
range of mass accretion rates 



a.d v 



/Cat for a 



Condition CvAheat ~ 1 at L 



- crit 
J u, corc 



To a factor of ~ 2 the 
However, the specific 
value of the ratio changes from about 0.7 to almost 1.2 
over the range of mass accretion rates - more than a 
50% increase. Clearly, this definition does not exactly 
correspond to the critical luminosity. 

The se cond definition o f adve ction and heating times 
is due to IThompson et al.l ([200511) , who defined the times 
locally as 

Cv = ¥ . ( 20b ) 
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Figure 13. Top: Mass in the gain layer M ga i n (lower set of lines) 
and the total mass between the neutrinosphere and the shock Aftot 
(upper set of lines) as a function of L v , core. Filled circles at the end 
of the lines denote position of L" 1 * ore . Colors denote variations in 
the mass-accretion rate going from \M\ = 0.01 Mq s~ 1 (black) to 
\M\ = 0.58 Mq s (yell ow) in uniform logarithmic steps. Bottom: 
Energy a la Janka (2001) in the gain layer i? ga i n (lower set of lines) 
and the energy between the neutrinosphere and the shock Etot 
(upper set of lines). The meaning of symbols and colors is the 
same as in the upper panel. 




..... J° 51 cr 'g s s l ] 
Figure 14. Ratio of adv ection and heating times as defined by 



Murphy & Burrows (2008). Each line shows the ratio as a function 
of core luminosity L„ t corc and fixed M. The color coding corre- 
sponds to different values of \M\ going lef to to right from 0.01 Mq 
s _1 (black) to 1.95 Mq s _1 (red). The filled circles at the end of 
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Figure 15. Ratio of advection and heating times as defined by 
I Thompson et akl (2005) for two values of M. Each line corresponds 
to a radial profile of the ratio for a given luminosity going from 
1 X 10 51 ergs s _1 (black) to the critical value (red). Filled circles 
at the end of the lines mark the value just inside of the shock. 

P/p 



f" — 
"-heat 



(21b) 



where rp = (dlnP/dr) 1 is the pressure scale- height at 
a given radius. In Figure H31 we plot the ratio ^dvAheat 
as a function of radius for two values M. For the smaller 
M, the ratio spikes to high values below the gain ra- 
dius, while it stays consistently below unity above it, even 
though the critical luminosity was reached. Contrary to 
this, the higher M calculation tolerates ^advAheat > 1 m 
the gain layer for a number of sub-critical luminosities. 

Neither definition of advection and heating times yields 
the desired exact correspondence with L"' 1 c t ore . We tried 
to adjust the definitions in many ways, but failed to pro- 
duce a condition that would be valid over a wide range 
of parameters (e.g., M, M, and r u ). 

4.4. Antesonic condition 

In Section [5] we showed that the antesonic condi- 
tion (cy lt ) 2 /Ue SC = 0.1875 is equivalent to the crit- 
ical condition for isothermal accretion with pressure- 
less free fall at the shock, and it can be generalized to 
max(c|/'yg SC ) w 0.19r for a polytropic EOS. Does a sim- 
ilar condition hold also in the more realistic calculations? 
In order to test this, we studied the maximum value 
of the ratio of the adiabatic sound speed cs to the es- 
cape velocity, max(cg/ui_), similarly to the polytropic 
case. We plot in Figure [16] how this value changes with 
core luminosity and mass accretion rate. We find that 
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Figure 16. The antesonic condition. Maximum value of the ratio of the adiabatic sound speed eg to the escape velocity v eBC as a fu nction 
of core (individual lines). Colors distinguish between different mass accretion rates and have the same meaning as in Figure 1141 Filled 



to be a local condition on the sound speed in the accretion flow, our analysis in Section [2] and in Figures [6] and [7] shows that the quantity 
max (c^ / Vg BC ) is a merely a scalar metric for solution space of Euler equations and thus it is a global condition. 



the value of this parameter is surprisingly constant at 
max(c|/fg SC ) ps 0.20 over almost three orders of magni- 
tude in M at L„, CO re = ^" corc- Furthermore, we took 
all 1350 critical luminosities used to construct Figure H] 
and calculated a histogram of values of max (cg/v^ sc ) at 
L " core- We find that the histogram can be very well fit 
with a Gaussian with a maximum at 0.193 and with a 
width of 0.009, that is only 5% of its value! While from 
Figure[T6]we observe in critical values a slight trend with 
M, we also see from the same Figure that the numerical 
noise contributes significantly to the total scatteiEB 

How does the value of max(c|/w^ sc ) depend on in- 
put physics? We determined max(c|/^e SC ) for data of 
Figure [5] and found that there is some variability. For 
constant heating (solid red line) the ratio is about 0.15, 
while for no cooling (solid blue line) and heating propor- 
tional to r~ 3 (red dashed line) the mean value is about 
0.26. For all other cases of different physics, the mean 
value ranges from 0.19 to 0.22, quite close to the re- 
sult from our fiducial calculation. Furthermore, we also 



11 Although critical luminosities can be determined very pre- 
cisely with our code, other parameters like critical shock radii are 
determined less precisely. This is because \drg/dL Vi corel = oo at 
the critical luminosity, as can be inferred from Figure \S\ Hence, 
a small uncertainty in L" 1 < J orc is dramatically amplified in rg and 
related quantities, as can be noted in Figures j 141 and [TBI 



find that increasing T, the fraction of free-fall velocity 
of the incoming matter (eq. from 0.25 to unity in- 

creases the value of the ratio from 0.19 to about 0.22. 
We found in Section [5] that for polytropic EOS the an- 
tesonic ratio is linearly proportional to adiabatic index T. 
With this information at hand, the close agreement be- 
tween max (c|/Vg SC ) measured in our fiducial calculation 
and the value in isothermal accretion can be regarded 
as a coincidence. The fiducial calculation is not exactly 
isothermal (T > 1), which compensates for lowering T to 
0.25. Even so, despite the fact that L™* ore (M) changes 
dramatically with the large changes to the input physics 
considered in Figure [5J max(c|/w^ sc ) is virtually con- 
stant ^0.2. 

Another question is how does the radius at which 
c| / v% sc peaks relate to other significant radii of the prob- 
lem like the gain and shock radii? In the isothermal 
calculation this critical radius coincides with the shock 
radius, because at this position the escape velocity is 
smallest. In the full problem this is clearly not the case 
as shown in the lower left panel of Figure 2] Instead 
we find that radius of maximum c| /v^ sc nearly coincides 
with the gain radius, but that they are offset by as much 
as 2% in both directions. For the case of no cooling, the 
maximum occurs at the neutrinosphere, and there is no 
gain radius. As the radius of maximum c|/Wg SC closely 
tracks the gain radius in the fiducial calculation, their 
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d ifference might n ot be noticeable at all in simulations. 

iBurrows et al.l (|1995| ) proposed a similar "coronal" 
condition for the explosion, where T of the substantial 
amount of matter needs to be higher than the local es- 
cape temperature. However, we have demonstrated that 
the transition from accretion to wind does not require 
sound speeds above the escape speeds. Instead, the 
sound speed vs. escape speed condition is only a man- 
ifestation of the inability to satisfy the shock jump con- 
ditions. Indeed, as Figured] shows, the Be rnoulli integral 
is neg ative at explosion as noted also by IBurrows et all 
((I995j) . 

5. ANALYTIC TOY MODEL 

In this Section we discuss an analytical toy model 
that retains the important properties of the full 
numerical steady-state sol ution (as characterized b; 
Burrows fc Goshvlll993l and lYam asaki & Yamadai 120051 
20061) . namely the existence of £™' ore and two branches 
of solutions with different shock radii. We are able to 
derive explicit analytic expressions for L c J^ orc and rs 



Lv.core), then rj < 1. In the presence of heating we still 



(eqs. IB7I & IB8I) . The details of the construction of the 
toy model are given in Appendix [B] Here we discuss 
the most crucial points of the model and the results we 
obtain. 

5.1. Construction of the toy model 

Our toy model is based on the assumption of hydro- 
static equilibrium, conservation of energy, and the fact 
that i"coro occurs when the shock jump conditions can- 
not be simultaneously matched to any solution of the 
Euler equations (Section [5]). We find that we can real- 
istically model i^core usm g the solution to simple al- 
gebraic equations. Specifically, we obtain a quadratic 
equation for rs, which gives L£ rl ' ore when the discrimi- 
nant vanishes. The cost of this simplicity is that we must 
introduce a number of approximations, albeit reasonable, 
and are forced to make our model somewhat internally 
inconsistent. Therefore, our toy model produces only 
qualitative agreement with the numerical solution. Even 
this approach, however, enables us to grasp the basic el- 
ements of the physics involved and to extend them to 
more complicated and more complete models (see Ap- 
pendix [B]). 

A primary simplification of the model is that we take 
the heating and cooling functions to depend only on ra- 
dius r: H(r) = aL^ cmc /r 2 and C(r) = b/r 4 . While 
H is a good approximation, C lacks the important self- 
regulatory feature of realistic cooling, which is propor- 
tional to T 6 , that an excess in the internal energy can be 
quickly radiated away. In other words, when no heating 
is present, the material cannot radiate away via neutrino 
cooling more energy than its total internal energy con- 
tent. To quantify this constraint on the cooling function 
C(r), we assume that a cold parcel of matter falls from 
large distance through the shock to the neutrinosphere 
at r v . When no heating is present, the parcel will radiate 
a fraction of its gravitational potential energy 



^r 2 p r . GM 
: — L r dr = n . 

M r v ■ 



(22) 



consider equation ([22]) to be valid, but with a different 
(and potentially higher) value of rj. A condition similar to 
equati on (|22|) was employed by Fernande z fe Thompson! 
(2009) to limit excessive cooling due to discreteness ef- 
fects. 

Shock radii as a function of L VtCOTe and L™* ore as a 

function of M are shown in Figures [T7] and [T51 respec- 
tively, at two levels of approximation. In both, the re- 
sults from the numerical solution to equations (IB1|) and 
(|B2I) are labelled "numerical" and the results from a truly 
analyt ic to y mo del a re labelled "toy." We find that equa- 
tions (|B7P and (|B8j) from the toy model, which provide 
explicit analytic solutions for rs and I/™* ore , reproduce 
the qualitative behavior of the more complete solution 
remarkably well. 

5.2. Relation of L c J% rc to Lf c 

Is it possible to reach L" 1 c t orc with just L^ cc l Within 
the framework of the toy model, the answer to the ques- 
tion of whether L"^ mf} is always higher than rjGM\M\jr v 
depends on the ratio of the term in the square brackets 
in equation (|B10|) (which is always positive and larger 
than unity) to the denominator, 2irr u p v a, an upper limit 
of the optical depth t v for our fixed density profile. For 
small rj we can thus write (see eq. IB8[) 



- crit 



r) GM\M\ 



(23) 



where rj is the cooling efficiency. If the sole source of 
neutrino flux is the accretion flow itself (ie. there is no 



where the exact factor in front of the right-hand side is 
comparable to, but always higher than unity. For pa- 
rameters given in Table [T] we get r„ w 2irr u p l/ a = 0.36 
and i"c 0rG > 7]GM\M \jr v . We confirmed this result for 
all < rj < 1 by inspecting graphs similar to Figure [TH1 
From equation ([23]) we see that if t„ ~ 1, then even 
if the whole region between the neutrinosphere and the 
accretion shock is immersed in a neutrino flux corre- 
sponding to the maximum integrated neutrino cooling 
r/GMM /r„, the critical luminosity cannot be attained. 
Furthermore, if rj > r„ then the critical luminosity is 
even higher than the maximum gravitational energy re- 
lease rate. In reality, the accretion luminosity available 
for absorption will be substantially lower than in the ideal 
scenario described, and hence sufficient core neutrino lu- 
minosity independent of the accretion flow is required to 
reach an explosion. One could rightfully argue that the 
optical depth is not constrained within our toy model and 
that at high t„ some effect might come to play that low- 
ers the critical luminosity below the accretion luminosity. 
However, as demonstrated by our numerical solution to 
the more detailed version of the problem described in 
Section 13.41 and depicted in Figure [TTJ L"' 1 ^^ is in fact 
much higher than the maximum possible accretion lumi- 
nosity. 

5.3. ^"c 0rc as a function of dimension 

lOhnishi et all ((2006) and llwakami et all (|200l ) inves- 
tigated the shock radii and stability of 2D and 3D accre- 
tion flows with heating and cool i ng bas ed on steady-state 
models of lYamasaki fc Yamadai (|2005f ) . They found that 
going from ID to 2D increases the shock radii even when 
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Figure 17. Shock radii rg as a functio n of L„, CO rc for the toy model described in Section [5] and App endix FBI The solid blue line shows 
the analytic solution given by equation |[B7]|, the dotted green line is a numerical solution to equation l|B6j l where the depende nce of w e sc 
on rg was included, and the dashed black line is the num erical solution without fixing the density profile. All have rj = 0.4 (eq. [22]). The 
dash-dotted blue line shows the solution to equation ||B7| | with rj = 0.28. Critical points are marked with filled circles, and the upper and 
lower solution branches are shown with thin and thick lines, respectively. The inset plot shows density profiles from the numerical solution 

n52 „„„.„ „-l up tQ L£rit Qre (dashed purple 



^v, core • 



line). Solid lines for higher shock radii (upper solution branch) are for corc decreased in steps of 10 51 ergs s 1 down to the minimum 
for the upper solution branch. The dotted orange line shows the density profile assumed in the analytic toy model, p/p v = (r/r„)~ 3 . 
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Figure 18. Comparison of Iv n * ore from the numerical solution 
to eqs. JB1 b and JB2 ft (dashed blac k lin e) to the analytic solution 
for the toy model provided in eq. (IB8 1> (solid blue line; only the 
solution with lower L£, rl * ore is shown). The dash-dot red line shows 
r]GM\M\lr v (rj = 0.4; compare with Fig. IT?) . 



the region is convectively stable, and that this effect is 
more pronounced for higher L UlCOTe . They also noticed 
that their L VtCOIe = 6.0 x 10 52 ergs s _1 simulation ex- 
ploded even though it was stable in spherical symme- 
try and stead y-state calcula t ions. The common opin- 
ion shared by lOhnishi et al.l (|2006[ ) and Iwakam l et al.l 
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Figure 19. Illustration of effects of modified heating and cooling 
on rg (solid lines) and £v n * Dre (filled circles). Decreasing cooling 
(blue) and increasing heating (red) both act to decrease i^ rl * orc , 
and to increase rg at fixed L Vt corc with respect to the fiducial 
calculation (black). However, lower cooling efficiency at i^ rl * orc 
increases rg, whereas higher heating decreases rg at i^ rl * ore - 

(2008) is that the SASI is responsible for the shock re- 
vival, essentially lowering L™* orc . However, we have 
shown in our toy model (Figure [17] and in more involved 
calculations, Figure [8]) that decreasing the cooling effi- 
ciency rj lowers L™* ore , as displayed in equation ([23| . In 
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order to illustrate this point more clearly, we plot in Fig- 
ure [17] shock radii for somewhat lower cooling efficiency 
77. We see that not only is £"c ore decreased, but also 
that the shock radii consistently increase over the whole 
range of allowed core luminosities and that at fixed lumi- 
nosity the increase of rs is both relatively and absolutely 
more prominent. For lower 77 the shock radius at £™' orc 
increases only modestly, as can be seen from Figure [T71 
To put t hi s on more quantitative grounds, 
IQhnishi et ahl C 2 006) found non-exploding dynam- 
ical ID models for neutrino luminosities up to 
ii/,core = 6.5 x 10 52 ergs s _1 , while the models in 



2D did not explode up to L Vy 



5.5 x 10 52 ergs s 1 . 



These numbers can be identified as -L™* oro in ID and 



2D, respectively. Interpreted using equation (|23p . these 
values of L^ OTe suggest that 77 decreased by about 15% 
in going from ID to 2D. However, from the physics of 
the toy model and our derivation of equations (|23[) and 
(|B7|) . we find that changes in 77 are always accompanied 



by changes in rs at L" 1 *^, such that the ratio of two 
critical shock radii rg rlt for two cooling efficiencies 771 
and 772 at their respective L"^ OTC is 



Km 



Km 



(24) 



in their 2D simulations, yet L^ rl c t oro was lower than in ID. 
Furthermore, while both decreased cooling and increased 
heating give lower L^ mc , higher entropy, and higher rs 

at fixed L UyCOTC , the latter gives lower rs at -L™* orc , which 
we consider incompatible with results of hydrodynami- 
cal simulations. Finally, while it is somewhat difficult to 
disentangle the effects of modified heating and cooling, 
we think that the cooling is responsible for the reduction 
in i^ coro because the temperature dependence of cool- 
ing offers greater prospects for modifying the structure 
of the flow, providing lower L^ mc , higher rs, and higher 
entropy. Looking at Figure [5J we see that with no cool- 
ing £™* ore can be more than order of magnitude smaller 

at low |M|. As the usefulness of our ID calculations 
and the toy model are quite limited for discussing multi- 
dimensional effects, our work can provide predictions and 
scaling relations of various quantities that will be useful 
in understanding the results of future simulations. 

6. DISCUSSION AND SUMMARY 

We studied spherically-symmetric accretion flows in 
order to understand the physics of the critical neutrino 
luminosity, L£, rl * ore , which separates steady-state accre- 
tion from explosion — a formulation of the neutrino 
mechanism of supernovae. Our results can be summa- 
rized as follows. 



This equation implies that a decrease in 77 of 15% should 
correspond to an increase in rs at L"'^ OIC by 8 — 9%, 
in good agreement w ith the larger shock radii found by 
IQhnishi et al.l ([2006D in their 2D simulations (see their 

Fig. 4). 

Similarly, iNordhaus et al.1 (|2010( ) found that in 3D 
L"'^ mc drops to about 60% of the ID value (see our 
Figure IT2j) . What modification of the physics would be 
necessary to explain such a decrease? As shown in Fig- 
ure 1191 we find within our fiducial model that reduc- 
ing cooling by a factor of two yields the desired drop of 
L"' 1 * orc . Simultaneously, the entropy in the gain layer at 

fixed L^, core increases by about 30% and rs at £™' ore 
increases by about 5 to 7%. We obtain similar changes 
in -£/£, rl * ore and entropy if we increase the net heating 
q by a factor of two in the gain layer, thus simulating 
in a crude way increased energy deposition as might be 
provided by convection. However, in this case rg rlt de- 
creases by about 8 to 10%, al t hough at fixed L Vy corc rs is 
larger. lYamasaki fc Yamadal (|2006|) performed a calcula- 
tion similar to ours, which included convection in a phe- 
nomenological way, that essentially flattens the entropy 
gradient, and they found that while L"^ oie decreased 
substantially, the shock radii also decreased. 

We note here, that 77 is the only free parameter 
within our toy model and therefore the effects of all un- 
known physics get projected into 77. Even so, the de- 
crease in the critical luminosities and the increase of 
shock radii as a function of dimension observed in sim- 
ulations (e.g. [Ohnishi et al.l 120061: llwakami et al.l 120081 : 
IMurphv fc Burrowsll2008t INordhaus et al.H2010D may oc- 
cur because higher dimensional flows have less efficient 
cooling, and not bec ause of SASI or any oth er instability. 
This is supported by INordhaus et all (|2010ft who did not 
see in their 3D simulations as vigorous SASI activity as 



1. We investigated isothermal accretion flows with 
fixed mass accretion rate M and variable sound speed 
ct- In Figure Q] we showed that there is a maximum, 
critical sound speed c™ 1 * which allows for a steady-state 
shock standoff accretion shock. For ct > c^? 1 * the flow 
cannot conserve mass and momentum simultaneously at 
the shock and it must move outward on a dynamical 
timescale to establish a wind solution. For the shock 
jump conditions relevant for the supernova problem we 
showed that there is no smooth transition to an outflow 
solution. Wc studied the ratio c^/v^ sc and we find that 
it is constant at c™' and stays always below the value 
for the sonic point. Thus, for the first time, we derive 
the "antesonic" condition for explosion: if exceeds 
~ 0.197; 2 SC (eq. [4]), there is no solution for a steady- 
state accretion flow with a standoff accretion shock, and 
explosion necessarily results. The physical realization 
of this transition is likely modified by time-dependent 
hydrodynamical instabilities. 

2. The mechanism of c™' in the isothermal model 
does not depend on any particular heating mechanism, 
and we explicitly showed in Section 13.31 that it is 
directly analogous to L"^ ole observed in more complete 
steady-state calculations. Specifically, in Figure [7] we 
investigated and classified the structure of accretion 
flows in the supernova problem as a function of the core 
neutrino luminosity L v . coro with respect to the position 
of the sonic point and the constraints imposed by the 
shock jump conditions. We found that the critical 
solution does not correspond to the sonic point, but it 
is always marginally close (Figure [6]). We showed that 
the maximum of the ratio c|/w 2 sc at L™* ore is close to 
~ 0.19 (Figure ITBl and stays constant to within 5% over 
a wide range of M, and masses and radii of the PNS. 



22 



Pejcha & Thompson 



Thus, our "antesonic" condition corresponds directly to 
£"core an d is t nus superior in monitoring the approach 
to L c J l l mc when compared to other heuristic conditions 
proposed in the literature (Section [4]). For example, 
the ratio of the advection to the heating time in the 
accretion flow is not constant along the critical curve 
(Figures QI]& [15]). Instead, the physics of L c J* mc (M), 
and thus the neutrino mechanism of supernovae itself, is 
identical to the physics of the antesonic condition. 

3. Ex tendi ng the previous works of [Burrows fc Goshvl 
(|1993D and lYamasaki fc Yamadal (|2005j, |2006|) we studied 
the thermodynamical profiles of the flows for L v _ corc 
ranging from to £"'c 0rc (Figure @|. We confirm 
that there are two solutions with different shock radii 
and different energies for a given L„. coro (Figure [5]), 
which merge at i"c 0re - We also found that the critical 
curve robustly exists for a wide range of microphysics, 
although the exact normalization and slope vary (Fig- 
ure [8]). Finally, we quantified the dependence of ££ rl * oro 

over a wide range of M, and mass and radius of the 
PNS, providing a useful power-law fit as a function 
of the key parameters of the problem (Figure [9] and 
eq. [IS]). Our results imply that larger PNS radius 
decreases L"^ OIC , potentially favoring stiff high-density 
nuclear equations of state (see Section 3.4; eq. [18]). 
However, faster contractions of the PNS in soft EOS 
models can be favorable for explosion by increasing 
neutrino luminosities and energies. 

4. We include in our calculations a simple approximation 
to gray neutrino radiation transport (Figure [4j. We 
found that the neutrino cooling of the flow above 
the neutrinosphere, the accretion neutrino luminosity, 
L^ cc , lowers L c J^ mc by a small but significant amount 
(Figure ITT]) . Thus, it is the core, not the accretion flow, 
which is responsible for the transition to the explosion. 
Specifically, we show by analysis of our numerical results 
and by calculations within our toy model that £™* oro 
is always higher than the maximum available accretion 
power GMM/r v (Figure [H and eq. [23]). Even though 
-k™corc is always >L^ CC , we found in Figure H"2l that the 
inclusion of heating by accretion luminosity decreases 
L"'* orc by an amount comparable to the reduction 
caused by going from ID to 2D, or from 2D to 3D. 

5. Our numerical calculations and toy model imply 
(but do not prove) that the reduction of L"^ OTe seen 
in recent 2D and 3D simulations likely arises from an 
overall reduction in the cooling efficiency allowed by 
the accretion flows in multi-dimensions (Section 15.31 and 
Figure [I"?]) . and not from from additional heating caused 
by convection or the presence of shock oscillations. In 
particular, while both a reduction in cooling and an 
increase in heating increase the entropy in the gain 
layer and decrease L^ mc , a reduction in cooling always 

occurs with an increase in ^s(i"c re)' while an increase 
in heating results in a decrease in rs(L" 1 ' orc ). Because 
the simulations imply that ?"s(£™* orc ) is somewhat 
larger in 2D and 3D with respect to ID, our equation 



(|2~4"|) . derived from the toy model, implies that a decrease 
in the cooling efficiency of the flow is the cause of the 
observed decrease in i^ rl c t ore . However, more work on 
this issue is clearly warranted. 
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APPENDIX 

EULER EQUATIONS, EQUATION OF STATE AND NEUTRINO PHYSICS 

Using thermodynamic expansion of e and P we express the radial derivatives of v, p and T by inverting equations ([5]- 
[7]). The radial derivatives are 
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Combinations of thermodynamic variables in these equations are defined as 
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Our form of thermodynamic equations is the same as the Newtonian version of equations of iThompson et all (|200lD . 
except that we include dependency on Y e and we do keep Z and K separate, because they are related only for a special 
case of equation of state. 

Our equation of stat e includes relativistic el ectrons and positrons with chemical potential and nonrclativistic free 
protons and neutrons (Qi an &: Wooslevl I1996D . In our neutrino physics we use prescriptions of heating, cooling, 
opacity and reaction rates for charged-current processes^ with neutrons and protons from IScheck et all ff2006h . We 
explicitly calculate the degeneracy parameter of electrons and positrons, but we assume that degeneracy parameter 



We also modified the heating and cooling functions to in- 



clude electron-positron annihilation and creation, which increased 
in the fiducial calculation by about 2 to 5%. 
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Figure 20. Left: Absolute values of individual components of the Bernoulli integral for models from Figure[4] kinetic energy v 2 /2 (solid 
black), internal energy e (dashed red), pressure P/p (dotted green) and gravitational potential energy GM/r (dash-dotted blue). Right: 
absolute values of the radial derivative of the profiles in the top panel: v(dv/dr) (solid black), de/dr (dashed red), (dP/dr)/p (long-dashed 
green), P/p 2 (dp/dr) (dotted orange), and GM/r 2 (dash-dotted blue). 



of (anti)neutrinos is zero. We also assume that the energies of neutrinos and antineutrinos do not change between r v 
and rg. 

While we do not include nuclear binding energy in our energy shock jump condition (eq. [Tl]), we can assess its 
effects on rg and L^ ole with the help of Figure [Jj Nuclear binding energy would appear as a positive term on the 
left-hand side of equation (fTT|) . thus the total energy & just downstream of the shock would be smaller at a given 
radius. This corresponds to shifting the green line down in the energy panels of Figure [7] This would lead to smaller 
shock radii, because the common intersection for the momentum and energy jump conditions would occur at smaller 
radius. Similarly, L™* orc would increase, because the common intersection could occur at higher L v cmc . For example, 
consider the energy panel in lower left part of Figure [Jj If we shift the energy shock jump condition (green line) down 
by a small amount, it will intersect the flow profiles (black lines) exactly at the position of the similar intersection 
in the momentum panel (vertical dotted lines). Thus, a solution with a shock is possible and L c J^ ma increases. By 
moving the green line down even more, the intersection would move to lower radii, thus decreasing rg. We verified 
that these qualitative assessments are correct by a sa mple numerical calculation. Finally, we point out that nuclear 
binding energy was included in similar calculations of lYamasaki fc Yam ada (20(11). 

B. CONSTRUCTION OF THE TOY MODEL 

It is well known (e.g. Uank a 2001) that the region between the accretion shock and the neutrinosphere is very close 
to hydrostatic equilibrium. Figure [20] shows individual components of the Bernoulli integral (equation [13]) along with 
their radial derivatives. We see that contributions of both the kinetic energy term v 2 /2 to the total energy budget and 
of the v(dv/dr) term to the momentum equilibrium are negligible. Thus, we can formulate the steady-state problem 
without the v(dv/dr) term in the momentum equilibrium equation and neglect the v 2 terms on the left-hand side of the 
shock jump conditions. In this setting, the flow is hydrostatic and velocity can be formally defined as v — M /(4irr 2 p). 
We also set dY e /dr = and dL„/dr = 0, which implies t„ <C 1. With these simplifications the equations ([5H7]) can be 
reformulated as 

PW _ PH ,_r^ (bi) 

J r 2 r 

s ,/ x GM GM f 1 47rr' 2 p(r',r s ) .. , T 

h(n)-h(r 2 ) = +/ Q q{r',L UtCore )dr', B2 

n r 2 J r2 M 

where h = e + P/p is the specific enthalpy of the flow and we explicitly state that the density profile p(r, rs) depends 
on the shock radius rg . 

As boundary conditions we choose the two shock jump conditions assuming that the matter up the shock is in free 
fall, specifically h(rs) = GM/rs and P(rs) — v csc M / (Airr^) , and we fix the in ner b ound ary density to p(r u ) — p v . To 
make use of the boundary conditions, we set r\ — r v and r 2 = rg in equations (|B1I) and ()B2|) . We choose the equation 
of state of relativistic particles which satisfies e = 3P/ p and hence h = 4P/p. As discussed in Section [5] we choose 
q = H — C with heating and cooling functions defined as 

^ = ^P, (B3) 

C = A. (B4) 
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While H is a good approximation, C lacks the important self-regulatory feature of realistic cooling, which is proportional 
to T 6 , that an excess in the internal energy can be quickly radiated away. In other words, when no heating is present, 
the material cannot radiate away via neutrino cooling more energy than its total internal energy content. To quantify 
this constraint on the cooling function C(r), we assume that a cold parcel of matter falls from large distance through 
the shock to the neutrinosphere at r v . When no heating is present, the parcel will radiate a fraction of its gravitational 
potential energy 

^C(r)dr = ^, (B5) 
M t v 

where r\ is the cooling efficiency. If the sole source of neutrino flux is the accretion flow itself (ie. there is no L„, corc ), 
then X] < 1. In the presence of heating we still consider equation (|2"2")l to be v alid, but with a different ( and p otentially 
higher) value of 77. A condition similar to equation (f22j) was employed by iFernandez fc Thompson! (|2009f ) to limit 
excessive cooling due to discreteness effects. 

We verified with our relaxation code that our choices of equation of state and heating and cooling along with an 
assumption of hydrostatic equilibrium do not fundamentally change the phenomenology of the problem. Namely, even 
in this highly simplified version of the problem the critical neutrino luminosity and the two branches of solutions are 
still present. 

We construct our toy model by evaluating the momentum balance (equation |B1| ) and enthalpy conservation (equa- 
tion |B2| ) between the neutrinosphere and the shock and then solving for rs- In order to have an analytically tractable 
result we fix the density profile to be p(r) = p v (r /r„)~ 3 . However, as the problem is fully determined and the 
density profile is an integral part of the solution, an ad hoc prescription of the density profile violates one of the 
boundary conditions. We choose to keep the boundary conditions on enthalpy and pressure, h(rg) = GM/rs and 

P( r s) = Mv esc /(Trrg), and to violate the related condition on density p(rs) = Mv esc / (irGMrs). This internal incon- 
sistency means that our toy model will produce only qualitative agreement with the corresponding numerical solution. 
However, our point is to illustrate basic ingredients and physical principles of steady-state accretion, not to provide 
any numeric values. 

Integrating equations (|B1HB2[) and applying the boundary conditions we get a polynomial equation for the shock 
radius 

2wrlp v { aL v , core b \ ( r s ^ 4 | ( Mv esc ^Pvfl a^core j frs \ GM | np v b _ ^ 



M V K r v M 



where we neglect the dependence of the free-fall velocity v csc on the shock radius rs to get an analytic solution. We 
will discuss how this approximation changes the toy model results further below. With this assumption, equation (IB6|) 
is solved for the shock radius yielding 




(B7) 



While it is not possible to readily provide obvious physical interpretation to the individual terms in equation (|B7I) , we 
can say that they are combination of heating and cooling terms that include coefficients a and 6, which measure the 
strength of heating and cooling, the gravitational potential energy term proportional to GM/r u , and the ram pressure 
term that is proportional to v csc M. 

Our solution for the square of the shock radius has two critical points. The first critical point occurs when the 
denominator vanishes, L U: corc = b/2ar 2 . In this case one of the solutions stays finite and the other goes from negative 
(and hence unphysical) to infinite. This corresponds to the minimal neutrino luminosity that allows two solutions for 
the shock radiusQ Note that this happens for a finite and positive value of L u>COTe . 

The second critical point occurs when the discriminant vanishes. At this point the two solution branches merge 
and there is no solution for parameter combinations that give a negative discriminant. Thus, the neutrino luminosity 
corresponding to this critical point can be readily identified as the critical luminosity L£, rl * oro . The equation for the 
critical curve can be obtained by setting the discriminant in equation (|B7I) to zero and solving the quadratic equation 
for the neutrino luminosity 



aL 



crit 
za corc 



GM\M\ 
ixp v rl 



2-K 2 p 2 u rl 



± 



-M 



M 2 ^ GM\M\ / 1 _ 7t p u b v csc b 

GM M r u (GM) 2 irp v r v GM' 



Kp v r 2 v 



(B8) 



One might be tempted to interpret 4nr 2 pq/ Mdr 



as 

an equivalent and general condition for the minimum luminosity 
that gives two solutions for the shock radius. However, this is 



a mere coincidence arising from our choice of density pr ofile and 
^uation of state that cancels several terms in equations flBH and 
when rg — ► 00. 
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Table 1 

Values of parameters used in the toy model. 



Parameter 


Value 


Description 


M 


1.2 Mq 


neutron star mass 


TV 


40 km 


neutrinosphere radius 


pv 


3 X 10 10 g cm" 3 


density at the neutrinosphere 


^GSC 


-7.31 X 10 9 cm s" 1 


velocity just outside the shock 


V 


0.4 


cooling efficiency 


a 


6.45 x 10" 19 cm 2 g" 1 


heating function coefficient 


b 


GM\M\ 


cooling function coefficient 



where the physically relevant critical luminosity is the smaller of the two solutions. To explain this result, let us keep 
for the moment only the first two terms on the right-hand side. Merging the terms with a and b to g crlt and realizing 
that M = 4irrlpvv(r„), we can write q CIit r u /v(r u ) ~ GM/r v , and c| ~ qr v /v (similar to the discussion in Section [3^4]) . 
Combining these expressions shows that (c? 1 ') 2 ~ GM/r v , which is basically equivalent to the critical condition for 
isothermal accretion as given by equation @. The additional terms on the right-hand side arise because within the 
framework of the toy model we were able to calculate the critical condition exactly. 

We note at this point that the analytic and relatively short form of equations (|B6I IB8|) is a result of a particular 
choice of the density profile and the heating and cooling functions. It is indeed possible to repeat the derivations with 
different power-law indices for heating and cooling and even with a density profile prescription that satisfies the outer 
boundary condition on density, which would make the toy model internally consistent. The final polynomial equation 
would be of higher order and generally without analytic solution at all. Even without an analytic solution it is possible 
to determine the two critical points using Descartes' rule of signs (appearance of the second solution branch for the 
shock radius) and setting the polynomial discriminant to vanish (maximum neutrino luminosity that gives steady-state 
solution). However, we decide not to discuss these more complicated toy models as the lengthy equations are not more 
illuminating than the simplest version we present here. 

We set b = r]GM\M\/irp v as discussed in Section [5] and rewrite equation (|B8[) as 
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where the physically relevant solution is the one that gives lower critical luminosity. 

With values of other parameters of the problem given in Table Q] we can test how our toy model compares with 
our numerical solution with the same input physics. Figure [IT] shows the shock radius rs as a function of L v . corc for 

M = —0.5Mq s _1 . We compare the analytic solution given by equation (|B7j) (solid line) to the numerical sol utio n 
of the simplified problem (dashed line). With the dotted line we also plot the numerical solution to equation (|B6|) . 
which differs from the analytic solution by accounting for the fact that v csc is a function of the shock radius rg. Both 
analytic and numerical solution to equation (|B6[) agree very well within each other and hence we prove the assumption 
of constant w csc is justified. The agreement between the numerical solution without fixed density profile and the 
analytic solution is, however, only qualitative. From the inset plot, which shows the density profiles of the numerical 
solution, we see that the density profiles are to first order power laws with the slope dependent on the shock radius. 
This explains why there is a mismatch of ~ 8 x 10 51 ergs s _1 in the critical luminosity between the numerical and 
analytic solution. There is also a slight curvature to the numerically obtained density profiles that results in these 
shock radii always being slightly higher than the analytic ones. Despite the obvious inadequacies we are able to very 
well reproduce the numerical shock radii with our analytic model. 
Figure [T5] shows the comparison of the analytic critical curve given by equation ()B9|) with the numerical solution to 

the simplified problem. We see a very good agreement both in the absolute value and in the slope for \M\ < 1 Mq s _1 . 
At higher \M\, we see an upturn in both the analytic and numerical critical curves. In fact, the analytic critical 
curve ends when the term inside the square root vanishes. We were not able to get a smoothly extending numerical 
critical curve above ~ 2.5 Mq s" 1 , which suggests that it is absent or at least substantially changed at very high mass- 
accretion rates. However, we did not invest igat e this issue further as we do not see anything as nearly dramatic in the 
realistic critical curves presented in Section [33] and it is thus an artifact of our simplistic assumptions in the toy model. 
Furthermore, the mass-accretion rates where the critical curve ends are fairly high and steady-state approximation 
ceases to be valid in this area of parameter space. 

To make a direct comparison with the numerical results of the full fiducial model presented in Section [3.4l we expand 
the square-root term in equation 
efficiency (77 < 1) to get 

^crit 



B9[) in the limit of low mass-accretion rate (\M\ < lM a s 1 ) and low cooling 
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which is correct up to second order. The term in the square brackets is always positive and larger than 1. The 
denominator can be interpreted as an upper limit on the optical depth t„, because t v J Airpadr — 2i:r u p l ,a for our 
fixed density profile. 

We conclude with an explanation of derivation of equation (|24j) . The discriminant in equation (|B7j) vanishes at 
i"core) which allows us to plug in the expression for ££, rl * ore from equation (|B10[) . Under the same assumptions 
that were used to derive equation (|B10[) . we can calculate that the critical shock radius rg rlt does not depend on the 
particular L£ n * ore , but it is proportional to r\~ x l 2 . 



